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The  Freeman  model  is  a biologically  plausible  dynamical  neural  network  that 
mimics  the  complicated  wave  patterns  collected  from  the  electroencephalogram 
(EEG).  Freeman  has  studied  the  model  to  explain  how  biological  systems  interact 
and  extract  information  from  the  environment.  The  Freeman  model  has  hundreds 
of  free  parameters  that  need  to  be  fine  tuned  either  manually  (through  trial-and- 
error),  adaptation  algorithms,  or  numerical  simulations  to  match  experimental 
data. 

This  work  investigates  the  Freeman  model  and  its  dynamical  building  blocks 
as  the  basis  for  novel  information  processing  systems  based  on  nonlinear  dynamics. 
We  considered  three  issues.  First,  we  performed  a theoretical  analysis  of  the 
reduced  KII  (RKII)  basic  building  block,  where  simpler  dynamics  are  observed. 

We  showed,  using  the  center  manifold  theorem,  how  the  RKII  dynamics  are 
controlled  by  an  Hopf  bifurcation  and  designed  the  parameters  analytically. 

Second,  nonlinear  dynamics  of  coupled  RKII  sets  were  studied  analytically  using 
the  synchronization  of  output  channels  to  define  computational  primitives  based  on 
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specific  dynamics  of  the  network.  The  design  of  coupling  coefficients  in  the  network 
for  obtaining  a desired  output  response  is  based  on  our  theoretical  analysis  of  the 
RKII  set  dynamical  behavior.  The  computational  power  of  the  RKII  network  is 
demonstrated  by  two  applications  that  address  logic  computation  and  associative 
memory.  Thirdly,  we  designed  analog  VLSI  chips  to  implement  this  continuous 
dynamical  system.  Analog  design  faces  the  challenges  of  power  consumption 
and  limited  chip  area.  To  build  a more  efficient  integrated  circuit,  the  nonlinear 
function  in  the  KII  architecture  is  eliminated  by  redesigning  the  summing  node.  In 
the  new  design  the  summation  block  serves  simultaneously  as  a current  adder  and 
nonlinear  activation  function  preserving  the  original  functionality.  Chip  tests  show 
the  good  performance  of  the  new  design. 
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CHAPTER  1 
INTRODUCTION 

1.1  Olfactory  System  Model  as  Information  Processor 

The  olfactory  system  is  one  of  the  oldest  parts  in  the  brain.  It  is  a well- 
developed  distributed  system  that  achieves  odor  identification  and  segmentation. 
The  olfactory  bulb  receives  inputs  from  specialized  odor  receptors  and  sends 
nonconvergent  signals  to  the  olfactory  cortex.  Olfactory  cortex,  in  cooperation 
with  other  parts  of  the  brain,  send  the  feedback  signals  to  the  olfactory  bulb  to 
accomplish  recognition  tasks.  Through  physiological  and  experimental  studies, 
many  researchers  have  built  computational  models  that  describe  the  architecture 
and  functions  of  olfactory  bulb.  In  this  thesis  we  will  focus  on  the  K-sets  model 
developed  through  the  seminal  work  of  Walter  J.  Freeman. 

The  realistic  computational  model  of  the  olfactory  system  proposed  by 
Freeman  [1]  describes  brain  function  as  a spatio-temporal  lattice  of  groups  of 
neurons  (neural  assemblies)  with  dense  interconnectivity.  It  is  a hierarchical  model 
that  quantifies  the  function  of  one  of  the  oldest  sensory  cortices,  where  there  is 
an  established  causal  relation  between  stimulus  and  response.  It  also  presents  the 
function  as  an  association  between  stimulus  and  stored  information,  in  line  with  the 
content  addressable  memory  (CAM)  framework  studied  in  artificial  neural  networks 
[2].  Freeman  uses  the  language  of  dynamics  to  model  neural  assemblies,  which 
seems  a natural  solution  because  of  the  known  spatio-temporal  characteristics  of 
brain  function  [3].  Although  we  believe  that  the  full  dynamical  description  of  the 
entire  model  (a  Kill  network)  is  beyond  our  present  analytical  ability,  one  may  still 
be  able  to  understand  the  dynamics  of  lower  level  structures  from  first  principles. 
One  advantage  of  a dynamic  framework  to  quantify  mesoscopic  interactions  is  the 
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possibility  of  creating  analog  VLSI  circuits  that  implement  similar  dynamics  [4,5]. 

In  this  respect,  dynamics  is  also  independent  of  hardware  (mimicking  the  well- 
known  hardware  independence  of  formal  systems).  However,  the  dynamical 
approach  to  information  processing  is  much  less  developed  than  the  statistical 
reasoning  used  in  pattern  recognition.  Only  recently  did  nonlinear  dynamics 
start  being  used  to  describe  computation  [6]  and  much  work  remains  to  achieve  a 
nonlinear  dynamical  theory  of  information  processing.  Hence  we  are  simultaneously 
developing  the  science  and  understanding  the  tool’s  capability;  a situation  that  is 
far  from  ideal.  The  challenge  is  particularly  important  in  the  case  of  Freeman’s 
model,  where  the  distributed  system  is  locally  stable  but  globally  unstable,  creating 
nonconvergent  (eventually  chaotic)  dynamics.  Nonconvergent  dynamics  are  very 
different  from  the  simple  dynamical  systems  with  point  attractors  studied  by 
Hopfield  [7],  because  they  have  positive  Lyapunov  exponents  [8].  Freeman’s 
computational  model  of  olfactory  system  has  already  been  applied  to  several 
information-processing  applications  [4,5,9,10].  Ultimately,  we  plan  to  use  Freeman’s 
model  as  a signal-to-symbol  translator,  quantify  its  performance  and  use  it  in 
analog  VLSI  circuits  for  low-power,  real-time  processing  in  intelligent  sensory 
processing  applications. 

Deterministic  systems  with  complex  dynamics  have  been  studied  in  the 
past  as  information  processors  and  feature  extractors.  In  addition  to  Freeman’s 
model,  other  approaches  using  oscillatory  or  chaotic  networks  [11-14]  have  shown 
great  computational  potential.  These  computational  systems  are  often  built 
upon  their  corresponding  dynamics  that  can  be  systematically  described  and 
controlled  by  energy  functions  or  convergence  properties.  To  accomplish  such 
a task,  we  should  construct  a well-defined  set  of  computational  primitives  to 
represent  the  state  space.  This  is  of  course  dependent  on  the  system’s  fundamental 
structures  and  its  specific  dynamics.  Different  dynamics  will  provide  a different  set 
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of  primitives.  Among  examples  of  information  processing  using  dynamical  systems, 
the  approach  is  usually  to  create  attractors  (either  convergent  or  nonconvergent) 
in  the  state  space.  The  case  of  dynamical  systems  with  point  attractors  (as  in  the 
Hopfield  network)  generate  local  minima  in  the  energy  function,  which  leads  to  the 
convergence  of  the  system  dynamics.  More  interestingly,  the  widely  used  recurrent 
networks  (RNN)  integrate  information  dynamics  with  system  dynamics,  which 
produces  better  representation  of  time  signal  and  more  efficient  projection  onto  the 
output  space  [15, 16].  Yet,  biological  systems  are  still  much  more  intriguing  in  both 
structure  and  observable  behaviors.  Brain  functions  deal  with  emerging  properties 
from  different  functional  layers,  where  nonconvergent  dynamics  are  the  carrier 
for  both  encoded  and  processed  information.  However,  the  case  of  nonconvergent 
dynamical  systems  (as  in  the  Freeman  model)  is  much  harder  to  quantify.  This  is 
especially  significant  for  the  Freeman  model,  where  detailed  resemblances  to  the 
real  biological  system  result  in  both  practical  implementation  issues  and  the  lack 
of  fundamental  understanding  of  its  dynamical  structures.  Existing  methods  used 
to  quantify  global  behaviors  of  the  network,  which  are  mostly  based  on  energy 
clustering,  are  still  far  from  satisfactory  to  unveil  the  advantages  brought  by  the 
complex  dynamics.  Fundamental  understanding  of  information  storage  (which 
depends  on  system  controllability  and  invariance)  and  information  retrieval  (which 
depends  on  the  understanding  of  projection  space)  both  remain  unsolved. 

One  possible  advantage  of  Freeman’s  system,  that  still  has  not  been  fully 
exploited  for  the  analysis,  is  the  uniformity  of  the  dynamics  of  its  individual 
components.  In  Freeman’s  model,  the  dynamical  behaviors  of  the  individual 
sub-systems,  such  as  fixed  point,  multiple  equilibria  and  limit  cycles,  are  very 
simple.  From  an  information-processing  viewpoint,  this  is  important  because  the 
global  dynamics  are  built  from  these  primitive  dynamics.  One  fundamental  step 
is  to  quantify  the  dynamics  of  the  individual  components  in  the  parameter  space. 
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Another  step  is  to  learn  the  rules  that  compose  the  complex  dynamics  of  the 
spatio-temporal  system  from  the  individual  dynamics.  This  is  challenging  because 
the  dimensionality  of  the  overall  system  can  be  very  high  to  lose  controllability 
and  invariance.  In  this  work,  we  used  a step-by-step  approach  with  the  goal 
of  unraveling  some  fundamental  mechanisms  that  naturally  emerge  due  to  the 
interconnectivity  of  the  components.  In  this  regard  we  used  synchronization 
analysis  among  the  processing  elements  (PE)  as  the  first  step  to  categorize  different 
regions  of  dynamics  in  the  parameter  space  to  elucidate  the  possible  composition  of 
global  dynamics.  In  fact,  global  behaviors  generated  by  interactions  among  lower- 
level  components  are  very  rich.  There  are  a wide  variety  of  collective  behaviors  in 
a coupled  oscillatory  system  including  synchronization,  partial  synchronization, 
desynchronization,  dephasing,  phase  trapping,  and  bursting  in  coupled  neural 
oscillators  through  diffusive  coupling.  These  are  all  possible  tools  to  refine  the  state 
space  of  our  system  to  provide  different  representations  of  embedded  signals.  If  all 
of  them  can  be  used,  computational  ability  and  efficiency  will  be  greatly  improved. 
However,  we  should  realize  and  emphasize  that  we  are  still  far  from  building 
computational  systems  upon  fully  explored  rich  dynamics,  where  parameter  design, 
controllability  and  system  invariance  need  to  be  thoroughly  investigated.  While  we 
could  not  provide  now  a full  interpretation  of  the  entire  system  dynamics  for  our 
model  , as  the  first  step,  we  tried  in  this  work  to  analytically  study  synchronization 
that  has  relative  simplicity  and  invariance  property  to  construct  an  initial  division 
of  dynamical  regimes  in  the  parameter  space.  By  controlling  one  of  the  emerging 
properties,  the  projection  space  can  be  built  so  that  in  these  complicated  and 
highly  interconnected  networks,  dynamics  collapse  to  a lower-dimensional  space  to 
facilitate  information  retrieval. 

Analytical  studies  have  the  added  advantage  of  decreasing  appreciably  the 
need  to  scan  parameters  that  may  still  be  needed  to  systematically  illustrate  the 
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global  dynamics.  To  achieve  the  ultimate  goal  of  designing  a desirable  compu- 
tational system  in  the  framework  of  a complex  dynamical  system,  it  is  of  great 
value  to  understand  these  basic  yet  rich  dynamical  behaviors  of  the  target  system, 
starting  from  the  most  basic  elements.  Moreover,  in  hardware  designs,  variations 
in  design  processes  and  fabrication,  as  well  as  noise,  may  degrade  the  performance 
of  the  system.  Theoretical  studies  help  recognize  the  range  of  the  system’s  working 
parameters  to  make  sure  that  variations  does  not  significantly  change  the  system’s 
qualitative  behavior. 


1.2  Structure  of  the  Freeman  Model 

Generally,  an  ./Vth  order  system  in  Freeman’s  model  is  defined  by 


1 
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where  a and  b are  experimental  constants  that  determine  the  cutoff  frequencies  of 
the  2nd-order  dynamics.  Wij  is  the  coupling  coefficient,  fj  models  a dispersive  delay 
and  Pi(t ) is  the  external  input.  Each  PE  in  Equation  1.1  models  the  independent 
dynamics  of  the  wave  density  for  the  action  dendrites  and  the  pulse  density  for  the 
parallel  action  of  axons.  Note  that  there  is  no  autofeedback  in  the  model.  Q(x)  is 
the  asymmetric  nonlinear  function  (at  the  output  stage)  in  each  PE.  It  describes 
the  pulse  to  wave  transformation  [1],  The  mathematical  model  and  properties 
of  Q(x)  are  discussed  later.  Freeman’s  model  is  a locally  stable  and  globally 
unstable  dynamical  system  in  a very  high  dimensional  space.  Although  it  has  great 
complexity,  the  model  is  built  from  a hierarchical  embedding  of  simpler  and  similar 
structures.  Based  on  the  seminal  work  of  Katchalsky  [1],  four  different  levels  (KO, 
KI,  KII  and  Kill)  are  included  in  the  model  defined  next  [1], 
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(a) 


(b) 


Figure  1-1:  Two  different  types  of  KO  sets,  (a)  Excitatory,  (b)  Inhibitory 


A KO  set  represents  about  103  to  108  neurons  that  receive  the  same  input  and 
have  the  same  sign  of  output.  Each  KO  set  consists  of  three  blocks  as  summation 
node,  2nd  order  linear  dynamics,  and  nonlinearity.  It  is  the  simplest  and  most 
basic  building  block  in  the  hierarchy.  All  higher-level  structures  are  made  of 
interconnected  KO  sets.  Spatial  inputs  to  a KO  set  are  weighted  and  summed.  The 
resulting  signal  is  passed  through  a linear  time-invariant  system  with  2nd-order 
dynamics.  The  output  of  the  linear  system  is  shaped  by  the  asymmetric  nonlinear 
function  Q(x).  Two  categories  of  KO  sets  (excitatory  and  inhibitory)  are  defined  by 
the  sign  of  the  nonlinear  function  (Figure  1-1).  There  is  no  coupling  among  any 
single  KO  set  when  forming  a KO  network. 

A KI  set  also  describes  a set  of  neurons  that  have  a common  input  source  and 
sign  of  output.  However,  there  are  dense  interconnections  within  the  KI  set.  KO 
sets  with  a common  sign  (either  excitatory  or  inhibitory)  are  connected  through 
forward  lateral  feedback  to  construct  each  KI  set  (Figure  1-2).  No  auto-feedback  is 
allowed  in  the  network. 

A KII  set  is  defined  to  realize  the  dense  functional  interconnections  between 
two  KI  sets.  Therefore,  there  are  three  possible  topologies,  which  are  formed 
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Figure  1-2:  KI  network  consists  of  interconnected  KOs  with  common  signs 


by  two  excitatory  sets,  two  inhibitory  sets  or  an  excitatory  and  inhibitory  pair, 
respectively.  The  last  one  is  specifically  denoted  as  a KII  set.  It  is  a coupled 
oscillator  consisting  of  two  KI  sets  (or  four  KO  sets)  of  different  signs  (Figure 
1-3).  Mitral  cells  Ml  and  M2  are  excitatory  cells,  while  granule  cells  Gl  and  G2 
are  inhibitory  cells.  Each  KII  set  has  fixed  coupling  coefficients  obtained  from 
biological  experiments.  A KII  set  is  the  basic  computational  element  in  Freeman’s 
olfactory  system.  The  measured  output  from  any  nonlinear  function  has  two  stable 
states  controlled  by  the  external  stimulus.  The  resting  state  occurs  when  external 
input  is  in  the  zero  state.  An  oscillation  occurs  when  external  input  is  present. 


8 


Therefore  a KII  set  is  an  input  controlled  oscillator.  The  KII  network  is  built  from 
KII  sets  that  are  interconnected  with  both  the  excitatory  cells  (denoted  by  Ml) 
and  inhibitory  cells  (denoted  by  Gl).  This  interconnected  structure  represents 
a key  stage  of  learning  and  memorizing  in  the  olfactory  system.  External  Input 
patterns  ( P(t ))  and  feedback  signals  through  Ml  cells  are  mapped  into  spatially 
distributed  outputs.  Excitatory  and  inhibitory  interconnections  enable  cooperative 
and  competitive  behaviors,  respectively.  The  KII  network  functions  as  an  encoder 
of  input  signals  or  as  an  auto-associative  memory  [1,4]. 

The  Kill  network  (Figure  1-4)  embodies  the  computational  model  of  the 
olfactory  system.  It  has  different  layers  representing  anatomical  regions  of  a 
mammalian  brain.  In  a Kill  network,  basic  KII  sets  and  a KII  network  are  tightly 
coupled  through  dispersive  connections  (mimicking  the  different  lengths  and 
thicknesses  of  nerve  bundles).  Since  the  intrinsic  oscillating  frequencies  of  each 
one  of  the  KII  sets  in  different  layers  are  incommensurate  among  themselves,  this 
network  of  coupled  oscillators  behaves  chaotically. 

1.3  Study  Objectives 

We  studied  Freeman  model  (an  oscillatory  network)  because  it  is  a computa- 
tion unit  that  mimics  biological  information  processing  with  sufficient  mathematical 
detail  to  enable  engineering  applications.  Freeman  wrote  extensively  about  this 
model,  but  kept  the  focus  of  explaining  biology,  not  to  help  build  man-made  infor- 
mation processing  artifacts.  Our  goal  is  to  theoretically  analyze  the  model  to  gain 
a basic  understanding  of  how  this  network  works;  how  to  control  it;  and  most  im- 
portantly, how  to  use  its  advantages  to  define  a system  of  computational  primitives. 
We  investigated  basic  levels  in  the  hierarchy,  emphasizing  dynamical  behaviors 
through  simulations  and  structural  analysis  that  explains  the  experimental  results. 
In  simulations,  our  main  task  was  to  identify  dominant  and  interested  dynamics 
from  different  behaviors  in  different  levels.  We  conducted  detailed  and  meaningful 
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Figure  1-4:  Kill  network  represents  the  computational  model  of  the  olfactory 
system 
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theoretical  analysis  for  lower-level  structures  from  the  K0  to  KII  networks.  Based 
on  the  theoretical  analysis  and  computational  primitive  defined  from  categorized 
dynamics,  potential  applications  were  investigated.  We  also  examined  hardware 
design  of  the  model  using  analog  VLSI  for  real-time  implementations. 

Chapter  2 begins  with  the  analysis  of  Hopf  bifurcation  in  the  reduced  KII 
(RKII)  set.  Hopf  bifurcation  characterizes  the  most  fundamental  dynamical 
behavior  of  the  basic  building  block  in  the  whole  hierarchy.  Analytical  solutions 
will  be  given  to  quantify  the  desirable  regions  of  control  parameters  in  the  model. 
The  results  will  then  be  used  to  provide  synchronization  analysis  of  higher  level 
structures  such  as  coupled  RKII  sets,  which  will  be  presented  in  Chapter  3. 
Chapter  4 discusses  the  computational  potential  of  RKII  networks  based  on 
synchronization.  Analog  VLSI  design  is  presented  in  Chapter  5.  New  architecture 
that  reduces  the  complexity  of  circuit  design  is  proposed.  Finally,  Chapter  6 
presents  conclusions  of  this  work. 


CHAPTER  2 

HOPF  BIFURCATION  IN  A REDUCED  KII(RKII)  SET 
A clear  understanding  of  the  dynamics  of  Freeman’s  model  is  needed  to 
control  the  network  and  to  define  its  characteristics,  so  that  we  have  insight  on 
how  and  why  to  use  the  model  in  different  applications.  Basic  observations  of 
dynamical  behaviors  of  the  olfactory  bulb  could  be  achieved  either  experimentally 
or  from  simulations  of  the  model.  However,  it  is  difficult  to  understand  the 
network’s  controllability  considering  the  massive  number  of  free  parameters  in  such 
a complicated  network.  With  a hierarchical  structure,  Freeman’s  model  is  built 
from  very  basic  components  such  as  a KII  set.  As  a general  approach,  the  study 
of  the  basic  component  helps  us  understand  both  local  and  global  behaviors  of 
the  network.  This  chapter  presents  a theoretical  solution  to  categorize  different 
behaviors  in  lower-level  structure  of  the  model. 

2.1  Dynamics  of  an  RKII  Set 


one  granule  cell  coupled 


Figure  2-1:  RKII  set  consists  of  one  mitral  cell  and 
through  Kmg(>  0)  and  Kgm(<  0) 
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An  RKII  set  is  a simplified  version  of  the  full  KII  set  [5]  shown  in  Figure 
1-3.  Unlike  a full  KII  set,  it  has  only  a single  excitatory  (mitral)  cell  and  a single 
inhibitory  (granule)  cell  (Figure  2-1).  The  M-cell  and  G-cell  are  coupled  by  two 
internal  coefficients  Kmg  and  Kgm.  P(t)  is  a time- varying  external  stimulus  but 
typically  it  only  has  binary  values  that  are  used  in  static  pattern  recognition 
[4,  5,  9, 10].  In  such  cases,  P(t)  is  assigned  either  an  off  state  (zero  input)  or  on 
state  (positive  input).  Although  an  RKII  is  structurally  different  from  a full  KII, 
the  network  dynamics  of  a large  assemble  of  RKII  sets  generate  similar  global 
dynamics  to  that  of  a full  KII  [5].  Moreover,  the  individual  dynamics  of  the 
RKII  set  are  qualitatively  similar  to  the  full  KII  set,  i.e.,  it  is  an  input  controlled 
oscillator.  For  our  analytical  studies,  the  RKII  has  the  great  advantage  of  fewer 
free  parameters,  which  simplifies  the  analysis.  In  addition,  the  most  significant 
property  of  the  olfactory  bulb  is  the  interaction  between  mitral  and  granule 
cells  (as  an  excitatory-inhibitory  network).  The  number  of  cells  is  normally  not 
balanced  and  their  connections  are  not  uniformly  distributed.  However,  starting 
with  the  simplest  yet  general  feedback  architecture,  higher  level  models  can 
always  be  achieved  by  well  defined  intracellular  connections.  From  a modeling 
perspective,  the  most  important  task  here  is  to  mimic  the  behavior  rather  than  to 
mathematically  match  the  target  biological  system.  Therefore,  we  focus  our  work 
on  the  study  of  an  RKII  set  and  hope  to  use  the  analysis  for  understanding  higher 
level  structures  in  the  system. 

An  RKII  set  is  described  by  a system  of  2nd-order  ordinary  differential 
equations  (ODEs)  as 

J_(^  + (a  + 6)^)+mW)  = KamQ{9)  + P 


(2.1) 
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where  m(t)  and  g(t)  are  state  variables  of  mitral  and  granule  cells,  and  Kmg  > 0 
and  Kgm  < 0 are  internal  couplings.  The  values  of  a and  b determine  the  cutoff 
frequencies  of  the  2nd-order  dynamics.  They  are  given  experimentally  as  a = 220  Hz 
and  b = 720  Hz  [1].  Q(x)  is  the  nonlinear  function  that  models  the  spatio-temporal 
integration  of  spikes  into  mesoscopic  waves  measured  in  the  olfactory  bulb  and 
serves  as  the  key  function  block  to  construct  this  population  model  [1].  Q(x)  is 
defined  by  Equation  2.2  as 

( Qm( l-e'W')  x > x0  (2.2a) 

Q(x)  = < 

{ -1  else,  (2.2b) 

where  x0  = ln(l  — Qmln(l  + 1/Qm)).  Qm  is  a free  parameter  that  determines 
the  ratio  between  positive  and  negative  saturation  values  of  Q(x).  In  our  study, 
for  simplicity,  we  used  Equation  2.2b  to  describe  both  positive  and  negative 
saturations  to  avoid  discontinuity.  This  variation  should  not  affect  our  results. 


Figure  2-2:  Plot  of  the  nonlinear  function  Q(x)  when  Qm  = 5 
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The  olfactory  bulb  is  a homogenous  excitatory  and  inhibitory  network  con- 
structed from  basic  building  blocks  (RKII  sets).  A typical  activity  of  an  excitatory- 
inhibitory  pair  is  the  birth  of  oscillations  induced  by  stimulus.  Actually,  in  such 
cases,  the  desired  dynamics  of  RKII  sets  are  expected  to  be  exactly  the  same  as 
those  of  the  full  KII  sets.  The  system  stays  at  its  equilibrium  (the  origin)  when 
external  stimulus  P(t ) is  zero  or  negative,  which  indicates  no  sensory  inputs  are 
received.  The  transition  into  oscillatory  states  (limit  cycles)  is  induced  by  a pos- 
itive and  large-enough  input.  A natural  approach  to  investigate  how  the  system 
undergoes  such  structural  changes  is  based  on  bifurcation  theory.  Many  well  stud- 
ied categories  of  bifurcations  led  to  the  birth  of  limit  cycles.  The  RKII  differs  from 
conventional  excitatory-inhibitory  structures  because  the  oscillator  itself  is  a 4th- 
order  system.  However,  next,  we  give  an  overview  of  dynamical  system  analysis  and 
show  that  the  transition  into  a limit  cycle  is  actually  through  a Hopf  bifurcation. 
This  reduces  our  current  model  due  to  the  fact  that  Hopf  bifurcation  resides  in  a 
2D  space  around  the  local  neighborhood  of  the  fixed  point,  yet  the  RKII  set  is  still 
biologically  plausible. 

2.2  Bifurcations  in  Nonlinear  Systems 

A complicated  nonlinear  system  is  usually  studied  by  linearizing  the  system 
equations  at  its  operating  points.  In  this  way,  its  local  behavior  can  be  described 
by  its  linear  counterpart.  Moreover,  the  dynamics  of  a linear  system  are  rather 
simple  and  are  dependent  upon  the  eigenvalues  of  the  linear  equations.  A linear 
time-invariant  system  is  defined  as 

^ = A x(t)  (2.3) 

where  x £ 5RA . A £ ?RNxN  is  a constant  matrix.  The  behavior  of  Equation  2.3  is 
determined  by  the  eigenvalues  of  A.  The  solution  x(t)  = x(0)eAt  has  the  following 
possible  dynamics  [17]: 
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• The  system  has  a stable  fixed  point  solution  if  the  real  parts  of  all  the 
eigenvalues  of  A are  negative. 

• The  system  is  unstable  if  at  least  one  of  the  eigenvalues  of  A has  a positive 
real  part. 

• If  every  eigenvalue  of  A has  a real  part  that  is  less  than  or  equal  to  zero,  and 
if  all  eigenvectors  corresponding  to  the  eigenvalues  with  zero  real  part  are 
independent,  than  the  system  is  stable;  otherwise  it  is  unstable. 

If  none  of  the  eigenvalues  of  A has  zero  real  part,  the  fixed  point  solution  of  Equa- 
tion 2.3  is  called  a hyperbolic  fixed  point.  Otherwise,  it  is  called  a nonhyperbolic 
fixed  point. 

The  analysis  of  a linear  system  is  rather  simple  and  straightforward  since 
all  possible  dynamics  are  clearly  determined  by  the  eigenvalues  and  eigenvectors 
of  A.  In  a nonlinear  system,  more  complicated  behaviors  may  exist  and  even  be 
unsolvable.  In  most  cases,  a nonlinear  system  is  linearized  around  its  equilibrium 
so  that  an  explicit  solution  and  qualitative  analysis  could  be  achieved  around  the 
neighborhood  of  the  fixed  point  [18].  Conditions  that  guarantee  a qualitatively 
similar  phase  portrait  between  a nonlinear  system  and  its  linearized  version  are 
described  by  the  Hartman-Grobman  Theorem  [19]. 

Theorem  1 (Hartman-Grobman  Theorem).  If  a nonlinear  system  x = f(x)  has 
a hyperbolic  fixed  point  at  x = xo  and  a Jacobian  A,  then  this  system  is  locally 
topologically  equivalent  to  the  linearized  system  (Equation  2.3). 

Based  on  this  theorem,  if  a nonlinear  system  (Equation  2.4)  has  no  purely 
imaginary  eigenvalues,  the  orbit  of  the  flow  generated  by  this  system  can  be 
mapped  to  that  of  its  linearization  as  defined  in  Equation  2.3  by  a homeomorphism 
(a  one-to-one  correspondence  of  points  in  the  two  flows). 

= f(x) 


dx 
d t 


(2.4) 
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Clearly,  stability  analysis  of  a nonlinear  system  can  be  greatly  simplified  while 
the  qualitative  property  obtained  remains  the  same.  However,  this  happens  only 
when  the  system  is  hyperbolic  and  not  at  a bifurcation  point  [19].  Bifurcation 
occurs  when  a system  is  structurally  different  (nonhyperbolic)  with  respect  to 
the  variation  of  its  parameter  set.  A parameter-dependent  system  may  present 
different  behaviors  in  the  phase  space  when  the  parameter  passes  through  the 
bifurcation  point  [19].  While  a nonlinear  system  can  be  reduced  to  its  linearized 
counterpart,  its  behavior  (or  trajectory)  around  the  bifurcation  point  cannot 
simply  be  described  by  a single  set  of  linear  systems.  In  fact,  various  bifurcations 
observable  in  biological  systems  are  usually  interesting  but  make  the  analysis 
difficult.  Some  basic  types  of  bifurcations  are  shown  in  Appendix  A.  The  question 
is,  when  a nonlinear  system  has  nonhyperbolic  fixed  points,  how  do  we  simplify  it? 
The  most  important  method  is  the  center  manifold  theory. 

2.3  Center  Manifold  Theorem 
Consider  a linear  system 

dx 

^ = Ax  (2.5) 

dt  ’ 

where  A is  a constant  matrix  and  x e UN . The  entire  space  9?^  can  be  described 

as  the  direct  sum  of  three  subspaces  determined  by  the  eigenvalues  of  A.  The  three 

subspaces  are  stable,  unstable  and  center  subspaces. 

• if  {ei,  ...,es}  are  the  generalized  eigenvectors  of  A whose  corresponding 
eigenvalues  all  have  negative  real  parts,  the  stable  subspace  is  defined  as 

Es  = span{ei, ...,  es}. 

• if  {es+i,  ...,es+u}  are  the  generalized  eigenvectors  of  A whose  corresponding 
eigenvalues  all  have  positive  real  parts,  the  unstable  subspace  is  defined  as 


Eu  = span{es+i, ...,  es+„}. 
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• if  {es+u+1,  ...,es+u+c}  are  the  generalized  eigenvectors  of  A whose  corre- 
sponding eigenvalues  are  all  purely  imaginary,  the  center  subspace  is  defined 
as 

E spon{es-|-u_|-i , es+u-|.c}. 

The  three  subspaces  are  all  invariant  manifolds,  which  means  that  the  solutions 
of  the  linear  system  with  initial  conditions  in  any  of  the  subspace  will  stay  in  it 
forever. 

Given  a system  without  any  solution  in  the  unstable  space  (that  is,  Eu  = 0),  it 
is  clear  that  any  orbit  or  solution  will  converge  to  the  center  manifold.  Depending 
on  the  dimension  of  the  center  manifold,  asymptotic  behavior  of  the  bifurcating 
system  under  analysis  can  be  approximated  by  the  reduced  system  locally  around 
the  fixed  point.  This  is  the  basic  idea  of  the  center  manifold  thorem  [17]. 

Theorem  2 (Center  Manifold  Theorem).  The  nonlinear  system  defined  by  2-4  can 
be  parameterized  by  states  in  the  center  subspace  (i.e.,  x G Ec)  as 

C = {x  + y|  y = 0(x)}.  (2.6) 

y = </>(x)  G Es  U Eu  is  called  the  center  manifold,  is  a smooth  function  and 
satisfies  <j>( 0)  = 0 and  0)  = 0. 

Furthermore,  the  center  manifold  is  locally  invariant  with  respect  to  the 
original  system,  i.e.,  every  solution  of  the  nonlinear  system  (Equation  2.4)  that 
originates  from  the  center  manifold  will  stay  in  the  center  manifold.  The  center 
manifold  can  also  be  attractive.  When  the  unstable  subspace  does  not  exist  (i.e., 

Eu  = 0j,  all  solutions  of  the  original  system  that  is  close  enough  to  the  center 
manifold  will  approach  the  center  manifold. 

Next,  we  characterize  the  dynamical  properties  of  the  RKII  set  and  use  the 
center  manifold  theorem  to  reduce  its  dimension  of  interesting  dynamics. 
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2.4  Hopf  Bifurcation  in  an  RKII  Set 

By  changing  of  variables,  Equation  2.1  can  be  described  as  a system  of  first- 
order  ODEs  by 


d m(t) 


= -abg  - (a  + b)gv  + abKmgQ(m), 


= —abm  — (a  + b)mv  + ab(KgmQ(g)  -I-  P) 


(2.7) 


where  m and  g represent  the  state  variables  of  mitral  and  granule  cells  respectively. 
mv  and  gv  are  the  lst-order  derivatives  of  m and  v.  The  internal  couplings  Kmg, 
Kgm,  and  input  P(t)  are  the  possible  bifurcation  parameters  in  this  system.  Before 
looking  into  the  the  actual  causes  and  properties  of  the  bifurcation,  we  first  discuss 
some  of  RKII’s  properties. 

2.4.1  Properties  of  the  Equilibrium 

The  equilibrium  point  can  not  be  solved  analytically,  but  we  will  investigate 
its  existence  and  uniqueness  with  the  help  of  numerical  simulations.  Let  us  assume 
m*  and  g*  are  the  the  fixed-point  solutions.  Based  on  Equation  2.7,  Equation 
2.8  can  be  used  to  solve  numerically  for  the  actual  value  of  the  equilibrium.  The 
equilibrium  should  be  at  the  intersections  of  the  two  nullclines  defined  by  the  two 
equations  in  Equation  2.8. 


m - KgmQ(g)  - P = 0 
9 KmgQijYl)  ~ 0 


(2.8) 


First,  we  discuss  the  the  existence  and  uniqueness  of  the  equilibrium  of  an  RKII. 
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Figure  2-3:  Equilibrium  is  determined  by  the  intersection  of  two  nullclines 

Property  1.  When  Kmg  > 0,  Kgm  < 0,  and  P > 0,  there  is  one  and  only  one  finite 
fixed-point  solution  in  an  RKII  set  defined  by  Equation  2. 7.  The  equilibrium  only 
exists  at  the  origin  or  in  the  first  quadrant,  i.e.,  m*  > 0 and  g*  > 0. 

Proof.  Recall  that  the  equilibrium  {m*,g*}  of  Equation  2.7  is  determined  by  the 
intersections  of  the  two  nullclines  m and  g defined  by  Equation  2.8.  The  system 
has  a fixed  point  at  the  origin  when  P = 0.  When  P > 0,  m and  g have  the 
same  sign  only  in  the  first  quadrant.  Hence,  the  intersection  is  always  in  the  first 
quadrant,  which  means  m*  > 0 and  g*  > 0.  When  P > 0,  it  is  straightforward  that 

Q'(x)  = exe~  > 0 


m'{g)  = KgmQ'(g)  < 0 
g'(m)  = KmgQ'(m ) > 0. 


(2.9) 
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Therefore,  m(g)  is  monotonically  decreasing  and  g(m)  is  monotonically  increasing. 
Therefore,  there  can  be  at  most  one  intersection  between  the  two  nullclines  m and 
g.  The  conclusions  can  also  be  verified  by  checking  the  intersection  between  the 
two  nullclines  in  figure  2-3.  □ 


Some  other  properties  of  the  fixed-point  solution  are  presented  below.  They 
will  be  of  interest  when  the  actual  values  of  the  internal  coefficients  need  to  be 
determined. 

Property  2.  Assume  that  m*  and  g*  all  have  finite  values,  when  P is  a constant, 
m*  is  a decreasing  function  of  both  Kmg  and  \Kgm\;  g*  is  an  increasing  function  of 
Kmg  and  a decreasing  function  of  \Kgm\. 


Proof.  According  to  property  1,  m*  > 0 and  g*  > 0,  so  we  have  Q(m*)  > 0 and 

ex  — 1 

Q(g*)  > 0.  Recall  that  Q'(x ) = exe  Q™  > 0.  The  proof  is  straightforward  based  on 
the  signs  of  the  lst-order  derivatives  computed  from  Equation  2.8. 


dm* 

* 

3 

* 

c? 

s 

1 

1 

dKmg 

1 + \KmgKgm\  Q'(m*)Q'(g*) 

dm* 

-Q{g*) 

d I Kgm  | 

1 + \KmgKgm\Q'(m*)Q'(g*) 

dg* 

Q(m*) 

dKmg 

1 -1-  \KmgKgm\ Q'{m*)Q'(g*) 

dg* 

-KmgQ(g*)Q'{m*) 

d | Kgm  | 

1 + \KmgKgm\ Q'(m*)Q'(g*) 

(2.10) 

□ 


2.4.2  Hopf  Bifurcation 

Intuitively,  the  behaviors  of  an  RKII  set  resemble  those  of  a Hopf  bifurcation 
during  the  onset  of  a periodic  orbit.  Specifically,  by  smoothly  adjusting  some  of 
the  free  parameters,  an  RKII  can  bifurcate  from  a stable  fixed  point  to  a stable 
limit  cycle.  This  matches  a supercritical  Hopf  bifurcation  (Appendix  A)  [17,  20]. 
Two  steps  will  be  taken  to  analyze  this  model.  First,  we  will  show  that  Hopf 
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bifurcation  exists.  There  are  three  free  parameters,  namely  Kmg,  Kgm  and  P. 

We  will  show  that  they  can  be  considered  together  as  a single  parameter  for  this 
particular  bifurcation.  Second,  properties  such  as  frequency  and  amplitude  of 
oscillation  of  the  periodic  solution  can  be  estimated  in  the  neighborhood  of  the 
unique  equilibrium  based  on  center  manifold  theorem.  Apparently,  we  expect  the 
center  manifold  to  have  a dimension  of  two. 

Conceptually,  the  Hopf  bifurcation  exists  in  a N dimensional  system  x = 

F(x,  /i),  where  x,  F € and  F 6 Cr,  r > 3,  when 

• Jacobian  of  F,  DF(0, //),  has  a pair  of  complex  conjugate  eigenvalues  A and  A 
such  that 

A (fj,)  = a(n)  + iu(fi)  (2.11) 


where  0)  = u>0  > 0,  a(g0)  = 0,  and  a'(g 0)  7^  0.  is  called  the  critical 
value. 

• The  remaining  N — 2 eigenvalues  all  have  strictly  negative  real  parts. 

To  examine  whether  a Hopf  bifurcation  exists  in  our  system,  we  first  compute 
the  eigenvalues  of  Equation  2.7.  The  corresponding  Jacobian  matrix  is 


A = 


0 

—ab 

0 

ClblPfYtgQ  {tyi  ) 


1 0 0 

■{a  + b)  abKgmQ'  (g*)  0 

0 0 1 

0 —ab  —(a  + b) 


(2.12) 


where  m*  and  g*  are  the  fixed-point  solutions  that  are  determined  by  Kmg , Kgm 
and  P.  Q'  (x)  is  the  derivative  of  Q(x ). 
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To  calculate  the  eigenvalues,  we  have 


|AI  — A|  = 


A 

ab 

0 


—abKmgQ'(m*) 


-1 

A + (g  + b) 
0 
0 


0 0 

-abKgmQ' (g*)  0 

A -1 

ab  A ~f~  (a  -f-  b) 


= 0. 


Thus,  we  obtain 


(2.13) 


[A(A  + a + b)  + ab}2  — (ab)2 KmgKgmQ  ( mn*)Q  ( g *)  = 0. 


(2.14) 


Solving  the  equation,  we  have 

A(A  + a + b)  — -ab  ± iabyJ\KmgKgm\  Q' (m*)Q' (g*).  (2.15) 


Therefore,  we  obtain  the  eigenvalues  as 


Al,2, 3, 4 — 


~(a  + b)±  J (a  - by  ± i4aby/\KmgKgrn\  Q' (m*)Q' (g*) 


(2.16) 


Note  that  Q'  (m*)Q'  (g*)  > 0. 


Clearly,  eigenvalues  of  A are  controlled  by  all  three  free  parameters  in  the  form 
of  g — \KmgKgm\  Q (m*)Q  (g*).  Without  loss  of  generality,  we  will  consider  these 
parameters  in  terms  of  a single  representation  by  g.  Due  to  the  fact  that  g,  is  not 
always  a monotonic  function  with  respect  to  any  of  the  free  parameters,  each  of 
them  may  lead  to  multiple  bifurcation  points.  Later  we  will  show  how  it  affects  the 
system  dynamics.  For  now,  we  only  consider  the  change  of  structural  stability  by  g 
in  a neighborhood  of  the  critical  value  giQ.  From  our  previous  discussion,  it  is  easy 
to  check  that  one  pair  of  the  system’s  complex  eigenvalues  always  have  negative 
real  part.  The  conditions  for  an  Hopf  bifurcation  are  all  satisfied  at  the  point 
where  a pair  of  purely  imaginary  eigenvalues  are  generated.  The  critical  values  are 
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Ho  = ( a + b)2 /ab  and  o>0  = y/ab.  In  a neighborhood  of  the  equilibrium,  the  frequency 
of  the  periodic  solution  can  be  approximated  by  the  system’s  time  constant  (1/a 
and  1 /b). 

To  determine  the  local  dynamics  of  this  bifurcation,  we  need  to  perform  center 
manifold  reduction  to  find  the  normal  form  in  the  neighborhood  of  the  equilibrium 
and  bifurcation  parameter  Ho  [17].  The  first  step  is  to  construct  the  real  Jordan 
canonical  form  of  A.  Transformation  matrix  T = [i?e{v!}  Im{vi}  i?e{v3}  /m{v3}], 
where  Vi  and  v3  are  the  eigenvectors  corresponding  to  the  pure  imaginary  eigen- 
value and  the  eigenvalue  with  negative  real  part,  respectively.  At  the  critical  value 
Ho,  we  have 


T_1AT  = 


(2,17) 


0 — y/ab  0 0 

y/ab  00  0 

0 0 — (a  + b)  — y/ab 

0 0 y/ab  — (a  + &) 

The  first  block  represents  the  invariant  center  manifold  while  the  second  block 
represents  the  invariant  stable  manifold.  Denote  y = T-1x.  The  new  system  y = 
G(y)  can  be  used  to  construct  the  normal  form  of  the  original  system.  By  solving 
the  center  manifold,  we  can  have  information  about  the  frequency,  amplitude  and 
stability  of  the  periodic  solution  in  the  neighborhood  of  the  equilibrium.  Analytical 
form  of  the  center  manifold  is  generally  unsolvable.  In  the  next  section,  we  will 
directly  analyze  the  bifurcations  as  a function  of  each  of  the  free  parameters 
to  obtain  accurate  solutions  that  can  control  the  dynamics  of  an  RKII.  Here, 
through  numerical  calculation  important  information  such  as  frequency,  amplitude 
and  stability  of  the  periodic  solutions  can  be  estimated  from  the  reduced  center 
manifold  [20].  We  should  notice  that  the  analysis  is  local.  Thus,  estimation  is 
accurate  only  in  a neighborhood  of  the  equilibrium.  From  simulations,  a reasonable 
estimation  can  be  obtained  in  a rather  large  range  beyond  the  critical  value. 
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Particularly,  the  periods  of  oscillations  in  a neighborhood  of  the  equilibrium  can  be 
estimated  using  ujq  as 

Ojr 

T = — [1  + k(a  - a0)  + o((a  - a0)2)],  (2-18) 

UJq 

where  cu0  = Vab  and  a0  = (a  + b)2/(ab).  The  values  of  k and  a need  to  be 
evaluated  after  normal  form  transformation  according  to  any  specific  set  of  system 
parameters  [20].  Figure  2-4  shows  an  example,  in  which  frequency  of  oscillation 


Figure  2-4:  Frequency  of  oscillation  as  a function  of  coupling  strength 

was  recorded  as  a function  of  /r.  Input  P — 0 and  Kmg  is  fixed  at  at  1 so  that 
H = Kgm-  The  bifurcation  point  is  at  (a  + b)2  jab  which  approximately  equals  to 
5.6.  Value  of  Kgm  is  increased  beyond  the  critical  value.  Estimated  frequency  from 
the  reduced  system  approximates  the  real  system  very  well  up  to  two  times  of  the 
coupling  strength  at  the  bifurcation  point. 
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2.5  Parameter  Design  in  an  RKII  Set 

Bifurcation  occurs  when  variations  of  parameters  result  in  structural  changes 
of  the  system.  In  other  words,  free  parameters  can  effectively  control  the  system 
behavior.  In  this  section,  by  studying  the  real  part  of  eigenvalues  of  an  RKII  set  as 
a function  of  coupling  coefficients  and  external  input,  we  can  find  ways  to  achieve 
desired  dynamics  of  RKII  sets  systematically. 

Recall  that  the  corresponding  Jacobian  matrix  is: 


A - 


0 10  0 
—ab  — (a  + b)  abKgmQ'  (g*)  0 

0 0 0 1 
abKmgQ'  (m*)  0 —ab  — (a  + b) 


(2.19) 


where  m*  and  g*  are  the  fixed-point  solutions  that  are  determined  by  Kmg,  Kgm 
and  P.  Q'  (x)  is  the  derivative  of  Q(x).  The  eigenvalues  are: 


Ai,: 


2,3,4 


-(o  + b)  ± \J (a  — b)2  ± i4aby/\KmgKgm\  Q'  (m\)Q'  (g\) 


(2.20) 


Let 

R = Re{)J ( a - b )2  ± i4ab^\KmgKgm\ Q' {m*)Q' (g*)} 

I = Im {\J (a  - b )2  ± i4abyj\ KmgKgm\ Q' (m*)Q' (g*)}  (2.21) 

It  is  clear  that  R > 0 and  (a  + b)  > 0.  Remember  that  system  behavior  switches 
when  Re{A}  passes  zero.  In  this  case,  — (a  + b)±R  determines  the  sign  of  Re{A},  so 
we  have  the  following  sufficient  and  necessary  conditions  for  the  two  distinguished 
states  of  the  equilibrium  that  we  are  interested  in. 

1.  An  RKII  set  (Equation  2.7)  has  an  unique  stable  equilibrium  iff 


R-max  < (a  + k)2 


(2.22) 
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2.  An  RKII  set  (Equation  2.7)  has  an  unique  unstable  equilibrium  iff 


Riax  > («  + bf 


(2.23) 


R = (a  + b)  is  the  bifurcation  point  of  the  system  that  represents  the  boundary  for 
the  purpose  of  controlling  the  behavior  of  an  RKII  set. 


Let  Rr  = (a  — b )2  and  Ri  = AabyJ\KmgKgm\  Q' (m*)Q' (<?*).  Note  that 


Rt 


Re2 { \jRr  4-  iRi} 


= sjR?r  + R2cos2(a/2) 
= 


1 + 


Rr 


yjm  + R\ 

yjRr  + R2  + Rr 


(2.24) 


By  solving  R^ax  = (a  + b)2,  we  obtain  the  bifurcation  point  in  terms  of  coupling 
coefficients  as 

(a  + b)2 


| Kmg  Kgm  | — 


1 


(2.25) 


Q'{m*)Q'(g*)  ab 

If  the  external  input  P is  fixed,  Equation  2.25  sets  the  boundary  between  stable 
fixed-point  and  limit  cycle  in  terms  of  Kmg  and  Kgm.  Generally,  the  right  side  of 
Equation  2.25  cannot  be  explicitly  solved.  However,  the  boundary  is  still  much 
easier  to  be  found  in  which  only  the  equilibrium  need  to  be  solved  but  the  whole 
parameter  space  doesn’t  need  to  be  scanned  and  checked.  In  the  case  of  zero  input 
(with  which  system  has  a fixed-point  at  the  origin),  Q'  (m*)Q'  (g*)  equals  to  1. 

The  values  of  a and  b are  predetermined,  so  we  have  an  explicit  upper  bound  as 
\KmgKgm\  = (a  + b)2 / ab  5.58,  which  agrees  with  the  critical  value  obtained  in  the 


previous  section. 
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In  general,  with  a fixed  external  input,  \KmgKgrn\  determines  the  dynamical 
behavior  of  an  RKII  set  by 

when  R2max  < (a  + b )2  (2.26a) 

when  Rmax  > («  + b)2.  (2.26b) 

Therefore,  for  the  purpose  of  controlling  the  dynamics  (i.e.,  fixed  point  solution 
with  zero  input  and  limit  cycle  with  positive  input)  of  the  system  defined  by 
Equation  2.7,  the  sufficient  and  necessary  conditions  on  the  coupling  coefficients  are 
• Condition  1:  if  and  only  if 

\KmgKgm\  < (2-27) 


\K-mgK gm\ 


1 


(a  + bf 


Kmg  Rgm  | 


Q'(m*)Q'(g*)  ab 
1 (a  + b )2 
Q'{m*)Q'(g*)  ab 


an  RKII  is  at  rest  when  there  is  no  input. 
• Condition  2:  if  and  only  if 


| R mg  Rgm  \ 


> 


1 (a  + bf 

Q'(m*)Q'(g*)  ab 


(2.28) 


an  RKII  oscillates  when  the  input  is  positive. 

Note  that  the  right  side  of  Equation  2.28  is  also  a function  of  Kmg  and  Kgm,  which 
is  nonlinearly  dependent  on  \KmgKgm\. 

Condition  2 also  sets  restrictions  on  the  value  of  Qm  in  Q'(x).  Q'(x ) is  defined 
by 

[ exe~^  X > x0  (2.29) 

Qf{x)  = l 

1 0 else  (2.30) 

where  x0  = ln(l  — Qm  ln(l  + 1 /Qm))-  Figure  2-5  shows  a plot  of  Q'(x).  Q'(x)  > 0 
and  has  a maximum  value  of  Qme~9^~  when  x = ln(Qm).  To  satisfy  condition 
2,  Q'(x)  must  be  greater  than  1 for  some  x.  Solving  the  maximum  values,  we  have 
Qm  > 1.  This  is  in  accordance  with  the  requirement  in  biological  models  that 
asymmetry  of  Q(x)  is  a necessary  source  of  instability  in  the  neural  network  [4], 
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Figure  2-5:  Derivative  of  Q(x)  when  Qm  = 5 


2.5.1  Determine  Proper  Values  of  Kmg  and  Kgm 

We  investigated  the  conditions  on  Kmg  and  Kgm.  In  this  section,  we  pro- 
vide the  specific  procedures  to  choose  suitable  Kmg  and  Kgm  that  satisfy  both 
conditions. 

(a  + b)2 


As  in  the  case  when  \KmgKgrn\ 


< 


ab 


-,  the  boundaries  of  Kmq  and  Kgm  is 


well  defined  by 


Kmg 

> 

0 

Kgm 

< 

0 

| Kmg  Kgm  | 

< 

(< a + b )2 

(2.31) 

The  second  case  is  more  complicated  because  the  2nd  condition  is  the  nonlinear 
function  of  Kmg  and  Kgm.  Although  we  could  not  calculate  an  analytical  solution 
of  the  region  where  valid  Kmg  and  Kgm  exist,  we  will  show  how  to  approximate 
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the  desirable  areas.  Again,  the  equilibrium  is  determined  by 


m*  - KgmQ(g*)  - P = 0 
g*  - KmgQ(m*)  = 0 


(2.32) 


where  m*  and  g*  are  the  fixed  point  solutions.  Solving  Equation  2.32  and  writing 
Q'(x)  in  terms  of  Q(x),  the  boundary  of  Condition  2,  \KmgKgm\  = 
could  also  be  defined  as 


where  m*(>  0)  and  g*(>  0)  are  also  functions  of  Kmg  and  Kgm.  It  is  hard  to  find 
an  explicit  solution  for  the  lower  boundary  of  desired  Kmg  and  Kgm,  especially 
in  the  form  of  \KmgKgrn\.  However,  we  can  still  find  ways  to  determine  useful 
properties  of  the  desired  regions. 

Property  3.  Given  m*  > 0 and  g*  > 0,  if 


(2.33) 


ab 


Q'{m*)Q'(g*)  = Q'  (m*)Q'(KmgQ(m*))  = 1, 


(2.34) 


then 


Q'(m*)Q'(KmgQ(m *))  > 1 


(2.35) 


for  all  0 < m < m* . 


Proof.  Q'(m*)Q'(KmgQ(m*))  has  similar  shape  as  Q'( x),  So  it  has  a lower  bound 
of  1 when  m decreases.  □ 
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Property  4.  Suppose  that  K^g  and  K*m  is  a set  of  values  that  satisfies  both 


gm 


conditions.  If  we  fix  the  value  of  Kmg  = K^g  and  increase  Kgm  until  \ K^K, 


(a  + bfi 

ab 


, then  \K  Kqm\  > 


'■mg 

1 


mg  gm  \ 


(a  + bfi 


m9  9m 1 Q'{m*)Q'(g*)  ab 


is  automatically  satisfied 


Proof.  At  the  point  where  I K*K* 


> 


1 


(a  + bfi 


, we  have 


' Q'(m*)Q'(g*)  ab 

Q'(m*)Q'(g*)  > 1.  With  increased  Kgm , both  m*  and  g*  are  decreased  (Prop- 
erties 1).  According  to  Properties  3,  Q' (m*)Q' (g*)  > 1 is  also  true  when  \K*K, 


(a  + b)2 

is  increased  to — . Thus, 

ab 


1 


\tc  K | = li±t£>- 
1 "*»  ab  Q'(m')Q'(g') 


(a  + bfi 


ab 


is  automatically  satisfied. 


mg  gm  | 

(2.36) 

□ 


Because  Equation  2.33  is  a continuous  function,  based  on  the  above  properties, 

/ j 2 

we  could  start  from  \KmgKgm\  = — , then  fix  the  value  of  Kmg  and  decrease 

the  value  of  \Kgm\  to  satisfy  both  Condition  1 and  Condition  2.  The  first  expo- 
nential term  in  Equation  2.33  will  be  monotonically  decreasing  while  the  other 
two  terms  will  be  monotonically  increasing  with  respect  to  a reduction  of  m*  and 
g*.  Two  cases  could  possibly  occur.  \KmgKgm\  could  be  dropping  until  it  reaches 


the  bifurcation  points  that  we  start  with.  Or,  it  will  reach  another  bifurcation 

1 ( a + b )2 


point  before  that.  Similarly,  starting  from  \KmgK, 


gm  | 


Q'{m*)Q'(g*)  ab 


and 


increasing  only  the  value  of  |A'9m|  will  also  give  us  a range  of  valid  values.  However, 
there  is  no  guarantee  that  these  two  areas  are  overlapped  so  that  starting  from 
any  valid  Kmg  and  Kgm,  all  values  with  larger  \Kgm\  will  satisfy  the  conditions. 
Nevertheless,  experiments  show  that  for  reasonable  and  normally  used  parameters 
(such  as  the  a and  b given  previously,  and  Qm—5),  we  can  use  the  above  procedures 
as  a rule  of  thumb  to  easily  determine  in  which  direction  we  should  change  the 
invalid  values  to  obtain  the  working  ones. 
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2.5.2  Dynamics  Controlled  by  an  External  Input 

The  state  of  an  RKII  is  controlled  by  an  external  stimulus  P.  Here,  P is  a 
time-invariant  input  with  two  states  ( off  and  on).  Off  state  (P  = 0)  is  guaranteed 
by  Condition  1 to  be  stable.  However,  not  every  positive  value  (on)  could  induce 
the  oscillatory  state  in  an  RKII  set.  We  will  see  in  the  following  that  P creates 
oscillatory  behavior  only  in  a bounded  region. 

Property  5.  In  the  system  defined  by  Equation  2.7,  the  fixed  point  solutions  m* 
and  g*  are  both  monotonically  increasing  functions  with  respect  to  the  external 
input  P. 


Proof.  The  lst-order  derivatives  of  m*  and  g*  with  respect  to  P are 

dm*  1 

dP  1 + \KmgKgm\ Q'(m*)Q'(g*)  > 

dg*  _ KmgQ'(m*) 

dP  1 + \KmgKgm\Q'{rn*)Q'(g*) 


(2.37) 


Recall  that  I\mg  > 0 and  Q'(x)  > 0 with  finite  x.  So  both  derivatives  in  Equation 
2.37  are  positive,  m*  and  g*  are  both  monotonically  increasing  functions  with 
respect  to  the  external  input  P.  □ 


Theorem  3.  Given  a set  of  fixed  coupling  coefficients  Kmg  and  Kgm,  the  values  of 
external  input  P that  enable  oscillatory  state  in  an  RKII  set  only  exist  in  a bounded 
region. 


Proof.  The  interval  of  x where  Q'(x)  > 1 (Figure  2-5)  is  limited.  The  conditions 
obtained  previously  require  Q' (m*  (P))Q' (g*  (P))  > 1.  Therefore,  m*  and  g* 
should  only  stay  in  a unique  finite  region.  According  to  Property  5,  P must  be 
bounded.  □ 
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0.15 


(b) 


o.i  - 


m(t)  VS.  g(t) 

* Initial  Value 


0.05 


0 


-0.05 


-0.1 


-0.25 


Figure  2-6:  Fixed  point  solution  in  an  RKII  set  when  P= 0.  (a)  Temporal  plot,  (b) 
Phase  plot 

2.6  Experiments 

2.6.1  Control  of  RKII  Dynamics  by  the  Coupling  Coefficients 

We  give  two  examples  in  which  the  external  input  is  fixed  to  either  zero  or  a 
positive  value.  In  the  case  of  zero  input,  an  exact  boundary  is  shown. 

When  P(t)  = 0,  Equation  2.26a  gives  a fixed  bifurcation  boundary  that  divides 
the  space  into  areas  of  oscillation  and  fixed  point  solutions.  As  given  in  [1],  a = 220 
Flz  and  b = 720  Hz.  So  when 


\KmgKgm\  < 5.5783, 


(2.38) 
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Figure  2-7:  Limit  cycle  solution  in  an  RKII  set  when  P— 0.  (a)  Temporal  plot,  (b) 
Phase  plot 

an  RKII  set  is  stable  and  has  a fixed  point  solution  of  (0,  0).  When 

\KmgKgm\  > 5.5783,  (2.39) 

an  RKII  set  oscillates.  Note  that  in  this  case,  the  actual  values  of  Kmg  and  Kgm 
do  not  change  the  qualitative  behavior  as  long  as  \KmgKgm\  is  unchanged.  In  the 
examples,  P = 0,  Kmg  is  fixed  to  one  and  Kgm  is  varied  to  set  different  values  of 
\KmgKgm\.  Initial  condition  of  m and  v is  set  at  (0.1,  0.1).  In  Figure  2-6,  when 
Kmg  = 1 and  Kgm  = —5.5  the  output  converges  to  a fixed  point  at  (0,  0)  after  a 
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Figure  2-8:  Fixed  point  solution  in  an  RKII  set  when  P=  1.  (a)  Temporal  plot,  (b) 
Phase  plot 

very  long  transient  time.  When  we  increase  \Kgm\  to  5.6,  the  RKII  set  converges  to 
a stable  limit  cycle  (Figure  2-7). 

In  the  second  example,  we  set  input  to  P = 1.  As  discussed  in  the  previous 
section,  there  is  no  explicit  expression  of  a lower  boundary  in  the  form  of  \KmgKgrn\ 
but  once  we  determine  one  set  of  valid  values,  all  other  values  with  a larger  \Kgm\ 
are  all  possible  choices  to  keep  the  system  oscillating  with  a positive  external  input. 
Again,  we  set  Kmg  to  1.  Figure  2-8  shows  the  transient  output  of  M-cell  and  G-cell 
with  an  testing  value  of  Kgm  = —4.  Apparently,  the  coupling  strength  is  not  large 
enough  to  sustain  an  oscillation.  In  this  case,  the  equilibrium  is  at  (0.1776,0.1906) 
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Figure  2-9:  Limit  cycle  solution  in  an  RKII  set  when  P= 1.  (a)  Temporal  plot,  (b) 
Phase  plot 


and 


| Kmg  Kgm  | — 4 < 


1 


(a  + 6) 5 


4.1852, 


(2.40) 


Q'(m*)Q'(g*)  ab 

which  is  less  than  4.  Trying  another  Kgm  = —5,  the  RKII  set  enters  the  oscillatory 
state  as  shown  in  Figure  2-9.  Here,  the  values  of  m*  and  g*  decrease  and  the 
equilibrium  is  at  (0.1502,0.1595).  Condition  2 is  satisfied  with 

1 ( a + b )2 


\KmgK, 


gm  | 


5 > 


Q'(m*)Q'(g*)  ab 


4.3762 


(2.41) 


After  scanning  different  values  of  Kgm , we  obtained  an  approximate  minimum  value 
of  4.237  as  the  lower  bound.  So,  in  this  case,  Kmg  = 1 and  —5.5783  < Kgm  < 
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Figure  2-10:  Fixed  point  solution  in  an  RKII  set  with  positive  input,  (a)  External 
input  P— 0.35.  (b)  External  input  P=26.1 

—4.237  are  all  possible  solutions  to  achieve  the  desired  dynamical  behavior  in  an 
RKII  set.  Note  that,  with  a relatively  small  \Kgm\,  the  amplitude  of  the  oscillation 
will  be  very  small  and  the  transient  time  to  reach  a stable  limit  cycle  is  rather 
large.  So,  normally,  \Kgm\  should  have  a relatively  large  value. 

2.6.2  Control  of  RKII  Dynamics  by  the  External  Input 

To  show  that  P only  exists  in  a bounded  region  to  achieve  desired  system 
dynamics,  we  fix  Kmg  = 1 and  \Kgm\  = 5 and  scan  different  values  of  P.  Approxi- 
mately, with  this  set  of  coupling  coefficients,  P must  be  kept  within  (0.42,  26.06)  to 
effectively  excite  the  RKII  set.  Figure  2-10  shows  the  transient  signal  of  the  system 
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Scanned  Boundaries:  P=1 


Figure  2-11:  Scanned  parameter  space  of  \KmgKgm 


with  two  different  P values  (0.35  and  26.1).  We  see  that,  not  all  positive  input  can 
force  an  RKII  into  the  oscillatory  state.  In  fact,  when  P = 0.35 

\Km9Kgm\  = 5 < Q'(m*)Q'(r)(  ^ = 5'0974’  (2’42) 

and  when  P — 26.1 

\KmgKgm\  = 5 < ^ = 5-2395’  (2‘43) 

2.7  Discussions 

The  Hopf  bifurcation  occuring  in  an  RKII  set  can  be  controlled  in  terms  of 
Kmg , Kgm  and  P.  The  condition  on  the  fixed-point  state  is  entirely  determined  by 

the  system’s  time  constants  and  can  be  solved  analytically.  Although  generally  the 

second  condition  on  input  induced  oscillation  does  not  have  a close  form  solution, 
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simple  rules  can  be  used  to  adjust  the  parameters.  As  long  as  the  first  condition  is 
still  satisfied,  the  coupling  strength  can  always  be  increased  to  achieve  oscillatory 
state  under  stimulation.  We  also  scanned  the  parameter  space  of  Kmg  and  Kgm 
to  give  a sense  of  the  shape  of  the  boundaries  (Figure  2-11).  The  simulation  is 
performed  with  P — 1 and  Qm  = 5.  We  see  in  the  plot  that,  the  lower  bound 
follows  the  properties  we  discussed  before.  If  an  arbitrarily  determined  coupling 
coefficients  could  not  make  the  RKII  set  controllable  as  we  expected,  we  need  to 
increase  the  value  of  either  Kmg  or  \Kgm\. 


CHAPTER  3 

SYNCHRONIZATION  OF  COUPLED  RKII  SETS 

A biological  system,  which  usually  presents  very  complicated  behaviors,  is 
often  a highly  interconnected  network  that  is  built  from  dynamically  similar  sub- 
systems. The  behaviors  of  such  sub-systems  can  be  very  simple  ones  such  as  stable 
fixed  point  and  multiple  equilibria  [21],  or  complicated  ones  such  as  oscillations  [22] 
and  even  chaos  [23].  An  interesting  question  is  how  the  complicated  behaviors 
of  the  overall  system  could  possibly  be  generated  by  the  interactions  among  its 
sub-systems.  This  leads  us  to  study  the  behaviors  of  coupled  systems.  There  are 
a wide  variety  of  collective  behaviors  in  a coupled  oscillatory  system  including 
synchronization,  partial  synchronization,  desynchronization,  dephasing,  phase 
trapping,  and  bursting  in  coupled  neural  oscillators  through  diffusive  coupling.  A 
broad  range  of  analysis  has  been  undertaken  from  physical  systems  to  biological 
systems  such  as  those  presented  in  [24-30]. 

An  important  aspect  of  this  research  work  is  to  understand  the  global  behavior 
of  Freeman’s  KII  network  so  that  we  can  achieve  an  efficient  and  robust  way  to 
use  the  network.  We  investigate  here  the  synchronization  of  the  basic  building 
block  in  the  computational  model  of  olfactory  system.  Basically,  Freeman’s  model 
presents  chaotic  behaviors  through  the  outputs  of  the  olfactory  bulb  layer  when 
it  is  either  in  the  basal  state  or  stimulated.  While  the  network  as  a whole  creates 
signals  that  are  similar  to  EEG  [31,32],  the  basic  building  blocks  (K0,  KI  and  KII 
sets)  are  either  single  or  coupled  nonlinear  dynamical  systems  that  have  simple 
dynamics  such  as  fixed  point  solution  and  limit  cycle.  Although  to  understand 
how  the  chaotic  behavior  is  generated  in  the  Kill  network  is  a very  challenging 
question,  by  using  synchronization  analysis,  we  expect  to  understand  one  of  the 
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global  behaviors  by  only  using  the  analysis  obtained  in  previous  chapters  to  avoid 
complicated  analysis  when  dealing  with  higher  level  structures.  It  can  be  a helpful 
first  step  to  start  understanding  higher  level  system  dynamics. 

Recall  that,  an  RKII  set  is  an  input  controlled  oscillator.  Given  the  external 
input  P,  a bifurcation  boundary  is  defined  by 


I KmgK, 


(' a + 6) 5 


grn  | 


Q'  (m*)Q'  (g*)  ab 

where  Q' (x)  is  the  derivative  of  Q(x)  and  ( m*,g *)  is  the  unique  equilibrium  of 
Equation  2.1.  Especially,  when  P = 0, 


(3.1) 


\Km,K,m\  = (3.2) 

In  this  case,  when  \KmgKgm\  < , Equation  2.1  is  stable  at  the  origin  while 

it  enters  oscillatory  state  when  \KmgKgm\  > . While  dynamical  behaviors  of 

individual  RKII  set  are  clearly  understood,  we  still  do  not  know  how  to  properly 
control  coupled  sets.  Analysis  on  synchronization  of  these  coupled  oscillatory 
elements  is  very  important  toward  our  goals  of  understanding  collective  behaviors 
of  higher  level  structures  in  the  Freeman  model.  Besides,  learning  and  memorizing 
in  the  olfactory  system  are  still  a challenging  problem  in  that  the  coding  of 
information  is  still  unclear.  The  OB  layer,  as  a coupled  RKII  network,  is  the 
key  element  for  information  encoding  and  computation.  It  has  been  used  as  a 
dynamical  associative  memory  just  like  the  other  oscillatory  networks  that  are 
being  used  in  learning  and  memorizing  [12-14].  The  analysis  of  synchronization  of 
coupled  RKII  sets  may  also  give  hint  on  possible  implementations  of  the  OB  layer 
(RKII  networks)  for  the  purpose  of  information  processing. 

Synchronization  of  coupled  chaotic  systems  has  been  widely  studied  [33,  34] 
due  to  many  challenging  applications  such  as  data  communications,  chemical 
reactions  [35-37]  and  collective  behaviors  in  biological  systems.  Pecora  and  Carroll 
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[33]  studied  the  synchronization  of  a subsystem  that  is  a replica  of  the  original 
system  where  the  driving  input  comes  from.  Later,  Ding  and  Ott  [38]  extended 
it  to  the  cases  that  the  synchronized  systems  are  not  necessarily  replicas  of  the 
original  system.  Brown  and  Rulkov  [39]  studied  the  sufficient  conditions  to  design 
proper  coupling  strength  of  two  chaotic  systems  that  guarantee  a synchronization 
between  them.  Although  the  analysis  are  for  coupled  chaotic  systems,  most  of  them 
also  apply  to  systems  with  simple  dynamics  such  as  a limit  cycle.  This  analysis  is 
of  particular  interest  in  our  study  of  neural  systems.  In  the  following  sections,  we 
discuss  how  to  use  Brown  and  Rulkov’s  method  to  determine  the  coupling  strength 
where  coupled  RKII  sets  synchronize. 

Section  3.1  presents  the  methodology  of  synchronization  analysis.  Analyt- 
ical solution  of  coupling  strength  to  achieve  synchronization  in  coupled  RKII 
sets  is  discussed  in  Section  3.2.  Simulation  results  are  presented  in  Section  3.3. 

The  derivation  of  the  analytical  solutions  of  the  sufficient  conditions  is  given  in 
Appendix  B. 

3.1  Sufficient  Condition  on  Synchronization  of  Coupled  Chaotic 

Systems 

The  synchronization  of  two  dynamical  systems  have  different  definitions  based 
on  different  purpose  and  applications.  Generalized  synchronization  (GS)  between 
a driving  system  and  a response  system  is  discussed  in  [40].  In  this  chapter,  we 
consider  the  case  that  follows  the  definition  of  identical  synchronization  (IS),  in 
which  the  transformation  from  driving  trajectory  to  response  trajectory  is  an 
identity  function.  Different  approaches  on  synchronization  analysis  of  coupled 
chaotic  systems  have  been  discussed  by  many  researchers.  Normally,  the  method  is 
based  on  Lyapunov  functions,  of  which  the  Lyapunov  exponents  are  difficult  to  be 
evaluated.  Brown  and  Rulkov  [39]  proposed  a sufficient  condition  on  linear  stability 
of  trajectories  in  the  synchronization  invariant  manifolds  of  two  identical  chaotic 
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systems.  In  their  method,  the  Lyapunov  exponents  do  not  need  to  be  explicitly 
calculated.  Also  based  on  rigorous  analysis,  they  gave  an  approximated  and  simple 
result  that  is  easy  to  be  calculated.  Although  solved  for  chaotic  systems,  the  results 
are  also  applicable  to  general  nonlinear  systems. 

Given  driving  and  response  systems  as 

dx  . v . . 

— = F(x;  t ) (3.3) 

and 

^ = F(y;()  + E(x-y)  (3.4) 

where  x,  y E 3?^,  x represents  the  driving  trajectory  and  y is  the  response  system. 
E is  a system  of  functions  describing  the  coupling  between  driving  and  response 
systems.  It  can  be  either  a linear  matrix  or  a set  of  nonlinear  functions.  Define  the 
deviation  of  system  y from  system  x as  w = y — x,  we  have 

r/w 

-=[DF(x;t)-DE(0)]w  (3.5) 

where  DF(i;  t ) is  the  Jacobian  matrix  of  F,  and  DE(0)  is  the  Jacobian  of  coupling 
matrix  E at  the  origin.  This  equation  describes  the  motion  transverse  to  the 
synchronization  manifold.  The  synchronization  manifold  is  linearly  stable  if  w (t) 
linearly  converges  to  zero  for  all  possible  driving  trajectories  x(t)  within  the  chaotic 
attractor  of  the  driving  system.  DF(x;f)  - DE(0)  can  be  further  decomposed  into 
a time  independent  part  A = (DF)  — DE(0)  and  an  explicit  time  dependent  part 
B = DF (x-,t)  - (DF). 

Assume  that  matrix  A has  eigenvalues  Ai,  A2, . . . , \n  and  are  ordered  by 
their  real  parts  as  3?{Ai}  > 5R{A2}  > . . . > Ji{A,v}.  Denote  the  matrix  of  the 
corresponding  eigenvectors  of  all  the  eigenvalues  as  T = [e1,e2,...,ew].  Brown  and 
Rulkov  [39]  solved  the  condition  for  linear  stability  of  the  invariant  trajectory  in 
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the  synchronization  manifold  as 


S[A1]><||T-1[B(x;f)T]||  >, 


(3.6) 


where  < • > is  the  operator  of  time  average.  In  addition  to  the  rigorous  result,  a 
quick  and  simple  rule  is  given  as 

»[Ai]  < 0,  (3.7) 

where  5R[Ai]  is  the  largest  real  part  of  the  eigenvalues  of  A.  That  is,  conditions  on 
synchronized  behaviors  can  be  studied  by  evaluating  the  eigenvalues  of  deviation 
system.  Next,  we  use  the  simple  rule  to  study  the  synchronized  behaviors  of 
coupled  RKII  sets. 

3.2  Analytical  Solutions  of  Couplings  for  Synchronized  RKII  Sets 


(a)  (b) 

Figure  3-1:  Two  representations  of  coupled  RKII  sets,  (a)  Linearly  coupled  RKII 
sets  (b)  Under  synchronization,  coupled  RKII  sets  could  also  be  described  as  an 
individual  RKII  set  with  auto-feedback 


By  the  original  definition,  coupling  among  different  KII  set  is  through  the 
same  nonlinear  function  Q(x)  as  defined  in  the  previous  chapter.  Linear  coupling 
between  oscillatory  components  in  olfactory  [41]  and  many  other  neural  system 
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models  has  also  been  widely  used.  While  different  couplings  all  address  the 
interconnection  among  different  set  of  neurons,  intuitively,  linear  coupling  has  much 
simpler  dynamics  and  is  easier  to  analyze.  Here,  we  shall  first  consider  the  case  of 
linear  coupling  through  coefficients  Kmm  and  Kgg  and  derive  an  analytic  solution. 

In  the  next  section,  we  discuss  the  effects  of  nonlinear  coupling  and  compare  the 
differences  and  similarities  between  the  two  coupling  scheme  in  terms  of  theoretical 
results  and  simulations. 

Assume  we  have  two  identical  systems  defined  by  Equation  2.1,  in  linear 
coupling  mode  (Fig.  3-1  (a)),  the  system  is  described  by  a system  of  ODE’s  as 


rn^ 


ml1'1 


abm(l)  - (a  + b)rriv)  + ab(KgmQ(g (1))  + P + Kmmrn (2)) 


g{1)  = gil) 


(3.8) 


gil)  = -abgW  - (a  + + ab(KmgQ(m^ ) + KgggW) 

where  Kmm  > 0 and  Kgg  < 0 are  the  excitatory  and  inhibitory  couplings 
respectively.  The  superscript  (1)  and  (2)  denote  the  two  systems  respectively.  In 
this  case,  coupling  matrix  DE(0)  is 


DE(0)  = 


abK„ 

0 

0 


0 

0 

0 


0 

0 

0 


0 

0 

0 


0 abKgg  0 


and 


A = (DF)-DE(O) 


0 


ofe(l  -(-  Amm) 

0 

abKmg(Q'  (m^)) 


-(a  + b)  abKgm(Q'(gW )) 


0 

0 


a6(l  + Kgg) 


0 

0 

1 

-(a  + b) 


(3.9) 


(3.10) 
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where  (Q'(x))  is  the  measure  (time  average)  of  the  derivative  of  Q(x).  The  outputs 
are  simple  limit  cycles  so  we  use  time  average  to  evaluate  them. 

To  find  the  sufficient  condition  that  guarantees  the  synchronization  between 
two  RKII  sets,  we  should  solve  for  the  real  part  of  the  Jacobian  of  A.  Detailed 
analysis  on  how  to  solve  for  the  conditions  will  be  presented  in  Appendix  B.  The 
results  are  given  here.  The  sufficient  condition  on  coupling  coefficients  Kmm  and 
Kgg  is  the  union  of  Equations  3.11-3.16. 

When 

Kmm  + \Kgg\  > 4 KB,  (3-11) 


either  of  the  following  three  inequalities  (Equations  3.12,  3.13  and  3.14)  needs  to  be 
satisfied. 


When 


(Kmm  - ?f)(\KM\  + ^f) 

A 

to 

Kmm  I K gg  | 

Km 

2 

(3.12) 

2 < Kmm  l-Rjgl  < 

Km 

2 

(1  + Kmm){\Kgg\  — 1)  < 

Kb 

(3.13) 

Kmm  1 Kgg  I 

Km 

2 

(1  + Kmm){\Kgg\  — 1) 

< Kb 

(Kmm  - ^f)(\K„\  + ^f) 

> Kb 

(3.14) 

0 < Kmm  + \Kgg\  < \J 4Kb,  (3.15) 
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Equation  3.16  needs  to  be  satisfied. 


-2  < K„ 


I K, 


99 1 


<!<„„„  4 I ft',, I)2  > 4AB  - 4A'.42  - 2KA'JKmm  - \K„\) 


(3.16) 


Values  of  Kai,  Ka2  and  Kb  are  calculated  as 

(a  - bf 


Kai  = 


ab 


Ka2 

Kb 


(a  + b)2 
ab 

KmlK,m(Q'(mm))(Q'(gW))  . 


(3,17) 


These  conditions  can  be  further  simplified  under  certain  assumptions.  Synchro- 
nized RKII  sets  could  have  different  kinds  of  behaviors.  For  our  purpose,  we  want 
these  coupled  sets  to  keep  the  same  dynamics  as  those  of  a single  RKII  set  when 
they  are  synchronized.  Under  identical  synchronization,  dynamics  of  two  coupled 
sets  can  be  described  as  a single  system  with  auto-feedback  (Figure  3-1  (b))  as 


rh  = mv 

rhv  = -ab(l  - Kmrn)m  - (a  + b)mv  + ab(KgmQ(g)  + P) 

9 9v 

gv  = -ab(  1 - Kgg)g  - (a  + b)gv  + abKmgQ(m).  (3.18) 

Kmrn  and  Kgg  interpret  couplings  from  the  other  system  as  auto-feedback  coef- 
ficients. Only  when  Kmm  < 1,  the  above  equation  defines  the  same  system  as  a 
single  RKII.  Thus,  we  limit  Kmm  to  be  less  than  1.  Moreover,  in  experiments,  the 
absolute  value  of  Kgg  is  normally  much  smaller  than  Kmm.  Thus  the  two  coupling 
coefficients  are  limited  to  the  unit  square.  Note  that  the  first  set  of  conditions 
require  Kmm  + \Kgg\  > 2 \JKB-  From  experiments,  KB  is  always  greater  than  1. 
Thus,  when  Kmm  < 1 and  \Kgg\  < 1,  the  sufficient  condition  on  synchronization  can 
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be  simplified  to 


{Kmm  + \Kgg\)2  > 4 KB  — 2KA2{2  + Kmm  — \Kgg\).  (3.19) 

3.3  Simulations 

We  consider  the  case  of  synchronization  when  an  uncoupled  RKII  set  is  in  the 
oscillatory  state.  Therefore,  if  a single  RKII  oscillates,  the  synchronized  systems 
should  also  stay  in  the  oscillatory  state.  In  order  to  keep  the  qualitative  behavior 
of  coupled  RKII  sets  the  same  as  the  behavior  of  an  individual  one  (oscillation  in 
this  case),  the  Jacobian  matrix  of  Equation  3.18  (i.e.,  DF(i*)  + DE(i’),  where  x* 
is  the  equilibrium)  should  have  at  least  one  eigenvalue  with  positive  real  part.  In 
this  case,  solutions  of  Kmm  and  Kgg  could  be  obtained  similarly  as  given  before  but 
with  a different  KB  as 


K* 


Kmg  Kgm  Q 


m 


)Q'(9 •) 


(3.20) 


Based  on  Equation  3.1,  an  oscillatory  RKII  set  must  satisfy 
KmgKgmQ'  (m*)Q'  (g*) 


= K“b>  E±fll  > 


ab 


(3,21) 


where  ( m*,g *)  are  the  equilibrium.  Thus,  similar  to  discussions  in  the  previous 
section,  we  obtain  a simplified  bifurcation  boundary  that  enables  oscillation  in 
synchronized  RKII  systems  as 


{Kmm  + \Kgg\)2  < 4 Kb  - 4 KA2  + 2 KA2(Kmm  - \Kgg\).  (3.22) 

We  denote  Equation  3.22  as  the  bifurcation  boundary.  Note  that  Equation  3.22 
is  only  valid  when  synchronization  is  achieved.  If  the  synchronization  condition  is 
satisfied  before  bifurcation  happens,  we  expect  the  coupled  sets  to  converge  to  a 
fixed  point.  The  bifurcation  boundary  is  more  accurate  because  it  uses  the  exact 
values  of  equilibrium  instead  of  estimations  of  limit  cycles  used  in  conditions  for 
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Figure  3-2:  Parameter  space  defined  by  Kmm  and  Kgg.  The  plot  shows  combined 
results  of  synchronization  and  bifurcation  boundaries  when  input  P = 0,  Kmg  = 1 
and  Kgm  = -6. 

synchronization.  The  two  boundaries  are  actually  two  parabolas  (Figure  3-2)  that 
divide  the  parameter  space  into  three  possible  regions,  where  two  coupled  RKII 
sets  behave  differently  as  synchronized  oscillation,  desynchronized  oscillation,  and 
fixed  point.  Figure  3-3  gives  an  example  of  a detailed  view  of  the  actual  boundaries 
inside  the  unit  square.  The  lower  bound  of  the  fixed-point  region  is  defined  by 
synchronization  conditions.  However,  we  used  the  time  average  of  limit  cycles  to 
estimate  Kb  but  in  the  fixed-point  region  the  trajectory  is  not  oscillatory  anymore. 
Therefore,  we  can  replace  KB  in  Equation  3.19  with  K*B  to  fine  tune  the  lower 
bound  of  the  fixed-point  region.  The  section  of  bifurcation  boundary  that  is  in 
the  desynchronized  area  is  removed.  There  are  three  other  parameters  (Kmg, 

Kgm  and  P ) that  could  change  the  position  and  shape  of  the  boundaries.  In  most 
cases,  these  three  parameters  are  predetermined  so  we  can  decide  the  boundaries 
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Figure  3-3:  Values  of  \Kgg\  and  Kmm  are  limited  by  1.  Three  regions  are  deter- 
mined where  the  coupled  system  shows  synchronized  oscillation,  desynchronized 
oscillation  or  fixed  point  solutions. 

accordingly.  In  the  case  when  their  values  need  to  be  changed,  it  is  still  convenient 
to  know  how  the  boundaries  move.  In  fact,  only  KB  (or  K*B)  will  be  affected  by 
the  three  parameters.  Increasing  or  decreasing  KB  (or  KB ) will  shift  the  parabolas 
along  the  position  vector  (1,-1)  (Figure  3-2). 

We  should  also  note  that  Equation  3.10  defines  a mutually  coupled  system. 
Driving  trajectory  is  not  determined  until  the  the  system  reaches  synchronization. 
This  affects  the  position  of  the  synchronization  boundary  because  the  value  of  KB 
changes  with  different  values  of  coupling  strengths  (hence,  the  synchronization 
boundary  may  not  be  strictly  a parabola  function).  If  KB  does  not  have  large 
variations,  it  can  be  simply  estimated  by  a single  oscillatory  RKII  set  without 
assuming  any  coupling  strength  in  advance.  To  achieve  a more  accurate  result, 
for  every  fixed  Kgg  ( Kmm ),  different  values  of  kmm  ( Kgg ) can  be  evaluated  by 
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computing  corresponding  values  of  KB  using  the  single  RKII  model  with  auto- 
feedback (Equation  3.18).  This  seems  to  lose  the  advantage  of  a simple  analytic 
solution.  However,  a single  RKII  set  usually  runs  much  faster  in  real  time  to 
reach  steady  state  for  a good  estimation  of  KB  than  two  coupled  RKII  sets  run  to 
achieve  synchronized  oscillation.  Single  RKII  set  is  also  simulated  more  efficiently 
because  it  has  fewer  dimensions.  It  should  also  be  pointed  out  that,  by  adaptively 
evaluating  the  coupling  coefficients,  only  the  synchronization  boundary  needs 
to  be  calculated  to  decide  the  sufficient  condition  on  synchronization.  This  is 
because  when  computing  the  time  averages,  we  assume  that  the  RKII  set  has  an 
oscillatory  state,  which  automatically  satisfies  the  bifurcation  boundary.  Moreover, 
when  external  input  is  zero,  the  other  two  boundaries  can  still  be  easily  calculated 
because  the  driving  trajectory  is  fixed  at  the  origin.  Apparently,  the  bifurcation 
boundary  should  be  always  below  or  overlapped  with  the  accurately  calculated 
synchronization  boundary. 

Next,  two  simulations  are  presented,  where  we  chose  Kmg  = 1,  Kgm  = -6 
and  P = 0 as  an  example.  First,  we  consider  two  cases  with  linear  coupling  where 
\Kgg\  = 0.8  and  \Kgg\  = 0.4,  respectively.  For  each  example,  we  show  that  only 
when  both  synchronization  and  bifurcation  conditions  are  satisfied,  the  coupled 
systems  behave  as  expected.  Second,  the  effect  of  nonlinear  coupling  is  discussed. 
Numerically  simulated  boundaries  and  calculated  boundaries  are  compared  for 
both  linear  and  nonlinear  coupling.  We  will  also  show  that,  in  addition  to  simple 
oscillations,  quasiperiodic  motion  is  also  observed  in  a small  region  around  the 
synchronization  boundary. 

3.3.1  Example  1:  Kgg  = 0.8 

When  increasing  Kmm  from  a very  small  value,  the  system  first  meets  the 
boundary  that  generates  synchronization  behaviors.  However,  in  the  synchro- 
nized system,  both  of  the  reduced  KII  sets  converge  to  a fixed  point  instead  of 
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Figure  3-4:  Simulations  of  coupled  RKII  sets  when  Kgg—~ 0.8.  (a)  Kmm- 0.6.  (b) 
Kmm= 0.96.  (c)  Kmm= 0.73.  (d)  Kmm= 0.74.  (e)  Kmm= 0.91.  (f)  Kmm  = 0.93 
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oscillating,  because  Kmm  has  not  passed  the  bifurcation  point  yet.  Theoretically, 
the  bifurcation  point  locates  at  Kmm  = 0.92  (Fig.  3-3).  Figure  3-4  shows  the 
steady-state  behavior  of  the  coupled  systems  at  Kmm  — 0.6,  0.96,  0.73,  0.74,  0.91 
and  0.93,  respectively.  When  Kgg  = -0.8  and  Kmm  - 0.6,  the  coupled  RKII  sets 
are  desynchronized.  When  Kmm  = 0.96,  the  coupled  RKII  sets  are  synchronized. 
When  Kmm  = 0.73,  desynchronization  occurs  below  the  synchronization  bound- 
ary. When  Kmm  = 0.74,  the  coupled  system  converges  to  a fixed  point  above  the 
synchronization  boundary.  When  Kmm  = 0.91,  the  coupled  system  converges  to  a 
fixed  point  below  the  bifurcation  boundary.  When  Kmm  = 0.93,  the  coupled  system 
synchronizes  and  oscillates  above  both  the  synchronization  and  bifurcation  bound- 
aries. The  results  are  as  expected.  We  see  that  theoretically  obtained  bifurcation 
boundary  predicts  very  well  experimental  results.  The  solution  serves  well  as  suf- 
ficient conditions  for  synchronized  dynamics.  Generally,  bifurcation  boundary  and 
lower  bound  of  the  fixed-point  region  are  more  accurate  because  they  are  calculated 
exactly. 

3.3.2  Example  2:  Kgg  = 0.4 

In  this  case,  the  bifurcation  condition  is  first  satisfied  when  changing  Kmm 
from  a small  value  (Figure  3-3).  However,  without  satisfying  the  synchronization 
conditions,  the  bifurcation  point  is  invalid.  In  this  case,  synchronized  behavior 
should  be  governed  solely  by  the  conditions  on  synchronization.  Figure  3-5 
shows  the  difference  between  the  two  systems  at  Kmm  = 0.33  and  Kmm  = 0.4. 
Synchronized  behavior  is  observed  when  Kmm  passes  0.37  while  theoretically  the 
sufficient  condition  is  about  0.42.  Again,  although  there  exist  errors  between 
calculated  and  actual  values  of  the  boundary,  the  analytical  solution  is  still  a 
sufficient  condition  on  synchronization  as  expected. 

Figure  3-6  illustrates  the  differences  between  analytical  solution  and  numerical 
simulation.  We  see  that  there  are  differences  between  the  two  synchronization 
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Figure  3 5:  Simulations  of  coupled  RKII  sets  when  Kgg  = —0.4.  The  coupled  sys- 
tem desynchronizes  and  synchronizes  at  (a)  Kmm  = 0.33.  and  (b)  Kmm  = 0.4 


54 


Figure  3-6:  Comparison  between  analytic  solution  and  numerical  simulation 

boundaries  but  the  analytical  solution  guarantees  synchronized  oscillations.  As 
discussed  before,  analytical  solution  gives  almost  the  same  values  of  the  other  two 
boundaries  as  those  in  numerical  simulation. 

3.3.3  Nonlinear  Coupling 

In  the  case  of  nonlinear  coupling,  Q(m)  and  Q(g)  are  coupled  through  Kmm 
and  Kgg.  Therefore,  weights  of  Kmm  and  Kgg  can  be  considered  as  time  varying 
according  to  a specific  trajectory.  To  use  the  simple  criterion,  similarly,  it  results 
in  estimating  time  averages  of  KmrnQ(m)  and  KggQ(g).  Therefore,  by  denoting 
K*mm  = Kmm{Q'(rn))  and  K*g  — Kgg(Q'(g )),  the  three  boundaries  (synchronization, 
bifurcation,  and  fixed-point)  can  be  rewritten,  respectively,  as 

(K-mm  + |ifs-s|)2  > 4 Kb  - 4 KA2  - 2 KA2(fCmm  - |jr;,|),  (3.23) 
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(Km m + NJ)2  < 4A-J  - 4^2  + 2KA2(K'mm  - \IC„\),  (3.24) 

and 

(K’mm  + \K'm\?  > 4 Ki  - 4 KA2  - 2KA2(K'mm  - (3.25) 

When  external  input  is  zero,  the  bifurcation  boundary  and  the  lower  bound  of 
the  fixed-point  region  should  be  exactly  the  same  as  those  in  the  case  of  linear 
coupling.  However,  simulations  show  that  the  simple  criterion  (Equation  3.7) 
fails  to  compute  the  synchronization  boundary  when  the  input  is  zero  and  all 
the  boundaries  when  the  input  is  positive  (where  in  nonlinearly  coupled  system, 
multiple  equilibrium  may  exist).  Table  3-1  compares  the  simulated  results  with 
analytic  solutions  when  Kmg  — 1,  Kgm  = —6  and  P — 0.  Data  shown  are 
the  boundaries  (synchronization  and  bifurcation  boundaries)  for  synchronized 
behaviors.  When  Kgg  = —0.7  (-0.8),  solution  of  the  lower  bound  of  the  fixed- 
point  region  for  both  linear  and  nonlinear  coupling  is  Kmm  = 0.681  (0.739)  which 
matches  the  numerical  simulation.  In  summary,  for  the  case  of  linear  coupling, 
boundary  for  synchronized  oscillations  can  be  calculated  using  the  synchronization 
condition,  where  with  a fixed  Kgg,  Kmm  can  be  evaluated  using  the  time  average 
estimated  from  a single  RKII  with  auto-feedback.  If  the  external  input  is  zero, 
part  of  the  synchronization  boundary  can  be  easily  calculated  using  the  bifurcation 
condition.  While  the  simple  criterion  generally  fails  in  the  case  of  nonlinear 
coupling,  it  still  applies  when  the  input  is  zero.  In  such  cases,  the  three  different 
dynamics  in  the  coupled  system  can  still  be  accurately  decided  in  part  of  the 
parameter  space  similar  to  the  case  of  linear  coupling. 

Table  3—1:  Comparison  between  analytic  solution  and  numerical  simulation 
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3.3.4  Quasiperiodic  Motion  in  Coupled  RKII  Sets 

Through  simulations,  it  is  found  that  the  coupled  system  shows  dynamic 
behaviors  consistent  with  the  parameter  partition  shown  in  Figure  3-2,  except  that 
quasiperiodic  motion  occurs  in  a narrow  region  below  the  actual  synchronization 
boundary.  From  simulations,  the  second  frequency  in  a linearly  coupled  system  is 
normally  slow.  So,  for  illustration  purpose,  a nonlinearly  coupled  system  is  used, 
where  Kmg=  3,  Kgm—  —3,  Kmm=  —0.4,  Kgg—  0.1885,  and  P=  0.  Figure  3-7  shows 
the  time-domain  signal,  the  phase  plot,  and  the  power  spectral  density,  respectively. 
In  fact,  the  transition  between  synchronized  and  desynchronized  oscillations  does 


Figure  3-7:  Example  of  quasiperiodic  motion  in  two  coupled  RKII  sets,  (a)  Tempo- 
ral signal,  (b)  Phase  plot,  (c)  Power  spectrum 

not  always  occur  with  a sharp  boundary.  Figure  3-8  shows  an  example  of  the 
region  where  quasiperiodic  motion  occurs  (between  the  actual  synchronization 
boundary  and  the  lower  bound  of  quasiperiodic  motion).  Note  that  this  behavior  is 
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only  observed  under  the  actual  synchronization  boundary.  Effort  was  made  through 
numerical  computations  to  find  out  whether  chaos  may  also  occur  for  linearly 
coupled  RKI1  sets  in  the  region  where  quasiperiodic  oscillations  are  observed.  Kgg 
and  Kmm  are  scanned  inside  the  region  where  quasiperiodic  motion  occurs  by 
fixing  Kgg  for  different  values  and  varying  Kmm  in  steps  of  10-4.  Power  spectrum  is 
examined  and  no  chaotic  motion  was  observed. 


Figure  3-8:  Simulated  regions  of  quasiperiodic  oscillations 


3.4  Discussions 

Many  biological  systems  are  modeled  as  interconnected  networks  of  neural 
oscillators.  Synchronized  oscillations  are  believed  to  play  a key  role  in  such  models, 
especially  for  the  purpose  of  information  processing.  In  applications  such  as 
associative  memory,  coupling  strengths  between  different  RKIIs  are  set  to  achieve 
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desired  temporal-spatio  response,  where  different  channels  are  clustered  based 
on  energy  and/or  phase  information.  Usually,  these  weights  are  set  empirically 
or  through  Hebbian  like  adaptive  learning.  In  adaptive  leanings,  proper  selection 
of  initial  condition  and  normalization  of  weights  are  normally  needed  to  avoid 
instability  in  training.  To  better  understand  the  dynamics  of  this  model  and  help 
for  network  design  in  applications,  we  studied  the  synchronization  of  two  coupled 
RKII  sets  based  on  a simple  criterion  to  gain  insight  on  how  to  properly  select  the 
coupling  coefficients. 

Linearly  coupled  RKIIs  can  be  analytically  solved  with  respect  to  kmm  and 
Kgg.  In  general,  with  fixed  Kgg,  the  sufficient  condition  for  synchronization  can 
be  solved  by  calculating  Kmm  based  on  the  oscillatory  trajectory  estimated  from 
a single  RKII  with  auto-feedback.  This  is  a simple  procedure  to  determine  the 
regions  of  different  dynamics.  It  can  effectively  avoid  simulating  the  coupled  system 
through  the  whole  parameter  space.  In  addition,  when  external  input  is  absent, 
part  of  the  synchronization  boundary  can  be  easily  obtained  due  to  the  stable 
fixed-point  at  the  origin.  The  criterion  does  not  apply  to  general  nonlinear  coupled 
systems.  However,  when  the  external  input  is  zero,  part  of  the  parameter  space  can 
be  solved  similar  to  that  in  the  linear  case.  More  sophisticated  method  is  needed  to 
fully  solve  the  problem  of  nonlinear  coupling.  However,  the  complexity  to  find  the 
boundaries  should  always  be  a concern  compared  with  the  effort  of  trial-and-error 
experiments.  Analytic  solutions  and  observations  were  validated  by  numerical 
simulations.  It  was  further  observed  that  near  the  synchronization  boundaries,  the 
dynamics  of  the  coupled  oscillators  may  become  quasiperiodic. 


CHAPTER  4 

COMPUTATION  BASED  ON  SYNCHRONIZATION 

In  this  chapter,  we  propose  to  use  synchronization  in  the  output  space  of  the 
Freeman  model  as  the  readout.  In  many  other  oscillatory  networks  [13, 14,  42], 
synchronization  among  groups  of  oscillators  has  been  used  in  various  applications 
such  as  image  segmentation,  audio  processing,  and  pattern  recognition.  More 
interestingly,  phase  coding  has  also  been  studied  for  Freeman’s  model,  which  shows 
advantages  over  conventional  pattern  recognition  techniques  based  on  amplitude 
[43].  In  this  chapter,  more  specifically,  we  will  not  only  use  synchronization  to 
construct  the  output  space  but  also  use  it  directly  in  designing  the  structure 
of  the  network  based  on  the  desired  functionalities  we  want  to  achieve.  In  this 
complicated  and  highly  interconnected  network,  in  response  to  input  stimulus,  the 
dynamics  collapses  to  a lower  dimensional  space,  and  different  groups  of  oscillators 
can  be  categorized  by  their  synchronized  behaviors. 

Based  on  the  discussions  in  previous  chapters  on  analyzing  dynamical  behav- 
iors of  the  basic  components  in  the  Freeman  model,  we  will  show  how  to  design 
critical  parameters  in  the  network  to  obtain  desired  output  responses,  especially 
synchronized  and  desynchronized  behaviors.  Two  applications  will  be  presented. 
The  first  application  provides  a different  view  of  computation  in  Freeman’s  model. 
With  fixed  coupling  coefficients,  a combination  of  external  inputs  will  lead  the  cou- 
pled oscillators  in  a reduced  KII  network  to  two  different  states  as  synchronization 
and  desynchronization.  It  turns  out  that  with  only  two  oscillators,  by  determining 
the  proper  values  of  the  coupling  coefficients,  we  can  build  functional  blocks  to 
implement  two-input  boolean  logic  gates  including  AND,  OR,  XOR,  NAND,  NOR, 
XNOR,  and  NOT.  In  the  second  application,  an  associative  memory  is  constructed 
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by  using  synchronization  among  output  channels  in  a reduced  KII  network  as  the 
readout.  The  network  combines  both  amplitude  and  phase  relations  among  RKIIs 
for  information  retrieval. 

4.1  Logic  Computation 

From  previous  discussions,  we  know  that  there  are  two  basic  behaviors  for  two 
coupled  RKII  sets.  If  we  denote  synchronization  as  state  1 and  desynchronization 
as  state  0,  the  coupled  system  is  a logic  function  when  defined  as 

/(I)  = S'(ro<1W2),0)  (4.1) 

where  I = {0, 1}  x {0, 1}  represents  the  input  space  and  S — {0, 1}  is  any 
function  that  determines  whether  or  not  the  outputs  of  the  two  systems  and 
m ^ are  synchronized.  0 is  the  threshold  used  in  S.  Note  that  I is  not  limited 
to  binary  values  of  0 and  1 . Basically,  any  two  values  that  represent  higher  and 
lower  level  inputs  can  form  the  input  space.  Function  S is  just  a symbolic  system 
where  one  can  assign  either  0 or  1 to  any  of  the  two  states.  This  means  that, 
if  / can  implement  one  set  of  the  basic  logic  computations  (for  example,  AND, 

OR,  XOR),  it  should  also  be  able  to  realize  their  complements  with  a different 
assignment  of  states  in  the  actual  measure  of  synchronization.  The  NOT  function 
can  also  be  implemented  with  many  options.  One  example  is  to  short  the  inputs 
of  NAND  gate  or  fix  one  of  its  inputs  to  1.  A two-input  logic  gate  has  a total  of 
16  boolean  functions.  Here,  we  only  consider  the  case  that  the  coupled  oscillators 
have  a symmetric  structure,  so  it  cannot  distinguish  between  the  inputs  (0,1)  and 
(1,0).  However,  this  may  be  solved  by  using  weighted  inputs.  The  weights  then 
become  the  parameters  that  need  to  be  designed  in  addition  to  Kmm  and  Kgg.  By 
considering  only  the  symmetric  structure,  without  combination  of  different  gates, 
the  coupled  oscillators  have  the  ability  to  implement  8 out  of  the  16  logic  functions 
that  include  all  the  boolean  primitives.  They  cover  most  of  the  universal  sets  (for 
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example,  {AND,  OR,  NOT}  and  {NAND})  that  can  express  any  logic  functions. 
Note  that  in  this  application,  analytic  solution  of  the  couplings  between  two  RKII 
sets  is  not  necessary.  Any  parameter  can  be  used  as  long  as  synchronized  and 
desynchronized  behaviors  can  be  achieved. 

Synchronization  can  be  measured  using  either  correlations  between  the  outputs 
of  the  two  sets  or  direct  detection  of  phase  difference.  Equation  4.2  calculates  the 
correlations  between  the  time  averagings  of  two  state  variables  . 

< m1m2  > — < mi  ><  m2  > 


C(mi,m2 ) 


(4.2) 


y/(<  m\  > — < mi  >2)(<  ml  > — < m2  >2) 

Notice  that  when  the  two  sets  are  synchronized,  the  phase  plot  is  a straight 
line.  When  they  are  desynchronized,  there  may  be  either  a limit  cycle  or  quasi- 
periodic  motions  in  the  phase  plane.  So  another  way  is  to  plot  the  bifurcation 
diagram,  so  that  multiple  period  could  be  detected  in  the  case  of  desyncronization. 

In  simulations,  we  use  nonlinearly  coupled  RKII  sets  and  fix  Kmg  = \Kgm\  = 3. 
Input  space  is  formed  by  binary  values  {0,1}.  Different  values  of  Kmm  and  Kgg 
are  searched  to  see  if  all  the  six  basic  logic  gates  could  be  realized.  We  relax  the 
limitations  on  the  two  parameters  and  allow  them  to  have  absolute  values  that 
are  greater  than  one.  Given  the  simple  structure  of  just  two  coupled  oscillators, 
the  results  are  very  promising.  Using  a synchronization  criterion  that  assigns 
state  1 to  synchronization  and  state  0 to  desynchronization,  the  coupled  RKII 
sets  can  implement  AND,  NOR  and  XNOR  functions.  That  is,  by  assigning  1 to 
desynchronization  and  0 to  synchronization,  we  have  the  complementary  functions, 
i.e.,  all  the  six  logic  functions.  Figures  4-1  and  4-2  show  the  implemented  XNOR 
function  in  the  time  domain  and  phase  plane  respectively.  Transient  behaviors  are 
discarded  to  guarantee  a steady-state  of  the  outputs.  The  outputs  of  the  nonlinear 
function,  which  have  saturations  at  -1  and  5,  are  plotted.  In  the  case  when  inputs 
are  [1  1]  and  [0  0],  the  two  reduced  KII  sets  are  perfectly  synchronized.  When 
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Figure  4-1:  Time  domain  response  when  implementing  coupled  RKII  sets  as  an 
XNOR  gate.  State  outputs  are  shaped  by  the  nonlinear  function.  Solid  line  shows 
Q(m^)  while  dashed  line  shows  Q(m^). 

inputs  are  [1  0]  and  [0  1],  the  two  oscillators  are  desynchronized  and  even  present 
chaotic  behaviors.  Note  that,  in  this  application,  zero  inputs  to  RKII  sets  still 
induce  limit  cycles.  Actually,  when  both  channels  are  synchronized,  there  are 
no  chaotic  behaviors  locally  or  globally  no  matter  what  the  inputs  are.  This  is 
different  from  the  conventional  Freeman  model  where  chaotic  behaviors  are  always 
present.  This  is  because  here  in  the  logic  computation  system,  the  output  space 
is  specifically  designed  to  be  either  synchronized  or  desyncrhonized  to  construct 
binary  functions.  And  as  discussed  previously,  when  two  channels  are  synchronized, 
we  expect  the  system  to  behave  as  a single  RKII  set  that  only  has  the  limit  cycle 
attractor.  Table  4-1  gives  the  parameter  values  in  two  coupled  RKIIs  to  build 
the  logic  gates.  In  the  second  column  of  the  table,  the  label  (0  or  1)  defined  for 
synchronized  behaviors  in  S is  indicated. 
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Figure  4-2:  Phase  plot  in  the  state  space  of  Mi  and  M2  when  implementing  cou- 
pled RKII  sets  as  an  XNOR  gate. 

Table  4-1:  Set  of  parameters  that  realize  the  complete  boolean  primitives 


Function 

Sync. 

Kmm 

K99 

Kfng 

Kgm 

AND  (NAND) 

1(0) 

1.75 

-2.0 

3 

-3 

NOR  (OR) 

1(0) 

1.2 

-0.4 

3 

-3 

XNOR  (XOR) 

1(0) 

1.5 

-0.4 

3 

-3 

One  important  question  is  how  to  extend  this  result  to  a network  of  oscillators. 
There  are  basically  two  structures.  The  first  one  is  to  implement  the  logic  gates 
just  as  in  digital  design.  That  is,  outputs  from  one  set  of  coupled  oscillators  will  be 
the  input  of  the  next  one.  And  combinations  of  different  gates  can  implement  more 
complicated  logic  functions.  The  other  way  is  to  fully  connect  N RKII  sets.  In  this 
network,  if  we  consider  the  synchronization  of  every  pair  out  of  the  N oscillators, 
there  will  be  a total  of  N(N  — l)/2  status  readouts.  Therefore,  the  size  of  the 
binary  output  space  could  be  as  large  as  2N<yN~1^2  (although  in  practice,  it  may  not 
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be  possible  to  design  a network  that  can  access  all  the  possible  outputs).  Although 
this  shows  another  aspect  of  computational  implementation  using  Freeman’s  model, 
we  should  also  notice  that  this  form  of  computation  still  falls  into  that  of  a normal 
system  due  to  the  much  simplified  binary  state  space. 

4.2  Associative  Memory  Using  an  RKII  Network 

An  auto-associative  memory  is  a content-addressable  memory  that  stores  a 
set  of  patterns  during  learning  and  recalls  the  learned  pattern  that  is  most  similar 
to  the  present  input.  There  are  basically  two  types  of  associative  memories:  static 
feed-forward  structures  and  the  recurrent  dynamic  structures  proposed  by  [11]. 
Hopfield  networks  have  a stable  point  attractor.  Here  we  will  be  showing  that  the 
reduced  KII  network  can  also  be  used  as  dynamic  associative  memories,  but  they 
differ  from  the  Hopfield  type  because  they  have  non-convergent  dynamics  (limit 
cycles  or  chaos). 

Let  us  consider  the  case  where  input  patterns  are  static  and  only  have  binary 
values  (positive  and  zero).  The  channels  (RKIIs)  that  receive  positive  (zero) 
stimulus  will  be  called  active  (inactive)  channels.  Ideally,  oscillatory  states  of  m’s 
that  are  receiving  the  same  inputs  should  be  grouped  together.  That  is,  we  expect 
zero  phase  lag  among  active  channels.  Under  a strong  constraint,  we  should  expect 
only  two  kinds  of  different  oscillations  that  are  corresponding  to  high  and  low 
input  values,  respectively.  However,  in  real  applications,  noise  in  the  background  or 
pattern  distortions  would  make  identical  synchronization  unrealistic.  Here,  we  will 
still  require  active  channels  to  be  synchronized  in  the  sense  that  phase  differences 
between  active  channels  only  deviate  from  one  another  in  a small  neighborhood 
around  zero.  Behaviors  of  inactive  channels  and  those  caused  by  errors  in  the  input 
space  will  be  relaxed.  Therefore,  the  only  criterion  for  recognizing  a pattern  is  to 
detect  the  active  group  that  has  strong  synchronization.  All  other  output  states 
that  are  not  strongly  correlated  to  this  group  will  be  classified  as  null-state. 
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Many  adaptive  methods  have  been  investigated  for  training  the  interconnection 
weights  (Kmm  and  Kgg)  in  Freeman’s  model.  Most  of  them  are  based  on  Hebbian 
learning  [4,  44,  45]  that  could  also  be  simplified  to  assign  binary  values,  strong  or 
weak,  to  coupling  coefficients  among  active  or  inactive  channels.  From  Chapter  3, 
we  know  exactly  how  to  control  the  behaviors  of  coupled  RKII  sets  by  different 
values  of  connections.  Thus,  coupling  strength  will  be  determined  targeting  specific 
desired  behaviors.  We  consider  four  scenarios  in  the  input  space  that  specify  the 
values  of  couplings  between  different  channels. 

To  construct  such  a network  as  an  associative  memory,  based  on  the  excitation 
coming  from  input  patterns,  the  following  cases  determine  the  coupling  strength 
among  different  channels: 

Active  channels  If  two  channels  are  to  be  excited  by  at  least  one  of  the  stored 
patterns,  when  applying  positive  stimulus,  we  expect  the  outputs  from  both 
channels  to  be  synchronized.  At  the  same  time,  when  receiving  incomplete 
patterns,  they  should  also  be  synchronized. 

Inactive  channels  If  two  channels  are  never  excited  by  any  of  the  stored  patterns, 
output  states  should  be  synchronized  with  weaker  oscillations  to  be  distin- 
guished from  the  active  channels.  At  the  same  time,  if  one  or  both  of  them 
are  excited,  they  should  be  able  to  desynchronize  from  one  another  so  that  all 
of  them  would  not  be  considered  as  active  channels. 

Active  and  inactive  channels  Without  the  overlap  between  different  channels, 
two  channels  that  belong  to  active  and  inactive  channels  respectively  within 
the  same  pattern  should  be  always  descynchronized. 

Overlapped  channels  If  two  patterns  share  the  same  channel,  connections 

between  this  mutual  channel  and  other  non-overlapped  channels  will  conflict 
for  different  patterns.  That  is,  for  one  pattern,  two  channels  are  both  active 
while  it  becomes  case  three  for  another.  In  such  situations,  intuitively,  we 
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will  reduce  the  inhibitory  connections  between  this  mutual  channel  and  other 
possible  active  channels.  Therefore,  the  activation  of  the  mutual  channel 
will  be  heavily  dependent  on  the  activations  of  other  active  channels  in  the 
incoming  pattern  but  not  affected  much  by  active  channels  from  another 
pattern  that  is  not  currently  presented  (inactive  for  the  incoming  pattern). 
Based  on  the  above  conceptual  design,  we  use  the  following  as  steps  to  set  coupling 
strength  for  an  RKII  network.  The  network  is  initialized  assuming  case  2.  For 
every  pattern,  the  coupling  strength  between  two  active  channels  is  set  according 
to  case  1.  Couplings  between  active  and  inactive  channels  are  set  according  to  case 
3 except  for  the  overlapped  channels,  which  will  be  set  according  to  case  4.  The 
four  cases  are  similar  to  the  process  of  Hebbian  learning.  However,  although  the 
couplings  are  designed  according  to  the  properties  of  input-output  correlations 
in  both  cases,  the  crucial  difference  lies  in  the  structure  of  the  readouts.  Outer- 
product  rule  or  adaptive  algorithms  such  as  Hebbian  learning  result  in  similarities 
among  mutual  channels  in  terms  of  energy.  However,  they  can  not  realize  our  goal 
of  guaranteeing  a solution  that  specifies  the  correlations  such  as  phase  information 
among  output  channels.  Moreover,  the  values  of  couplings  are  calculated  directly 
in  our  case  instead  of  through  learning  where  issues  such  as  convergence  rate  and 
stability  may  be  a concern. 

A 64-channel  RKII  network  is  implemented  as  an  associative  memory  that  is 
trained  to  memorize  and  recall  patterns  of  digits  “0” , “1” , and  “2”  in  the  size  of 
8 x 8.  To  differentiate  different  groups  of  synchronized  channels  and  their  levels 
of  activities,  we  use  the  following  quantitative  measure.  First,  cross  correlations 
Cy’s  between  the  ith  and  jth  channels  are  calculated.  Higher  correlated  channels 
are  grouped  together.  Denote  different  groups  as  Gi,  i < N,  where  N is  the  total 
number  of  channels.  The  evaluation  value  for  each  group  is  simply  described  by 
the  total  amount  of  input  received  by  all  channels  in  the  group  (Oj  = YlkeG  Ab 
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where  Ak  is  the  input  value  received  by  the  kth  channel).  The  same  evaluation 
value  is  then  assigned  to  all  channels  in  the  same  group.  Note  that  the  network 
is  designed  such  that  only  stored  patterns  will  be  recalled  as  synchronized  output 
channels.  Noise  comes  with  higher  input  amplitude  instead  of  zero  but  all  the 
output  channels  related  to  noise  or  false  input  will  mostly  not  be  synchronized  and 
not  be  grouped  together.  Of  course,  we  assume  that  noise  or  false  input  has  an 
amplitude  smaller  than  or  at  most  comparable  with  ideal  input  values.  For  most 
cases,  this  is  a simple  but  very  effective  criterion. 

Individual  RKII  set  is  configured  with  Kmg  = 1,  Kgm  = —6,  and  positive 
input  P = 3.  For  two  coupled  RKII  sets.  Case  1 couplings  are  set  as  Kmm  = 0.2 
and  Kgg  = —0.4.  Case  2 couplings  are  Kmm  = 0.2  and  Kgg  = —0.1.  Case  3 
couplings  are  Kmm  = 0.1  and  Kgg  = —0.8.  Case  4 couplings  are  Kmm  = 0.2  and 
Kgg  = —0.3.  All  couplings  are  further  normalized  according  to  number  of  channels 
in  the  network.  0.8  is  set  to  be  the  threshold  of  correlation  coefficients  between 
two  channels.  Two  channels  with  correlation  value  larger  than  the  threshold 
are  assigned  to  the  same  group.  In  the  presence  of  clean  patterns,  the  readout 
successfully  recalls  the  patterns  stored  in  the  network  (Figure  4-3(a)-(f)).  In  the 
time  domain,  we  see  that  all  channels  that  receive  positive  stimuli  are  synchronized 
and  have  much  higher  amplitude.  In  Figure  4-3(b),  desired  active  channels  are 
almost  identically  synchronized  (small  phase  differences  exist  between  the  active 
channels  for  “0”  and  those  overlapped  with  “1”)  and  have  larger  amplitude  of 
oscillation.  Desired  inactive  channels  have  small  amplitude  of  oscillation  and  are 
desynchronized  with  active  channels.  Using  the  criterion  given  in  the  previous 
section,  pattern  “0”  can  easily  be  recovered  to  the  stored  pattern  as  shown  in 
4-3(a).  In  Figure  4-3(d),  the  phase  differences  among  desired  active  channels  are 
slightly  larger  than  those  presented  in  the  case  of  pattern  “0”.  However,  using 
the  given  threshold,  pattern  “1”  can  be  perfectly  recovered.  In  these  cases,  a 
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Figure  4-3:  RKII  network  as  an  associative  memory,  (a)  Stored  pattern  “0”.  (b) 
Network  response  to  input  pattern  “0”.  (c)  Stored  pattern  “1”.  (d)  Network  re- 
sponse to  input  pattern  “1” . (e)  Stored  pattern  “2” . (f)  Network  response  to  input 
pattern  “2” 
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simple  energy  readout  would  also  work  perfectly.  However,  when  input  distortion 
(incompleteness  and  noise)  is  presented,  energy  or  distance  based  readout  may  fail. 


Figure  4-4:  Noise  is  added  to  an  inactive  channel  and  an  active  channel  in  “0” 

In  Figure  4-4,  positive  inputs  ( P = 3)  as  noise  are  fed  into  a desired  inactive 
channel  and  an  active  channel  for  “0”  only.  Figure  4-5  (a)  shows  the  total  amount 
of  input  received  in  each  synchronized  group  (Ot  = J^keGi  as  the  readout. 
Although  noise  increases  the  received  input  values  for  undesired  channels,  those 
channels  tend  to  be  grouped  separately  so  the  readout  value  will  be  kept  low. 

When  more  channels  are  distorted  by  noise,  as  long  as  they  are  desynchronized,  the 
final  readout  will  always  be  much  lower  than  the  value  of  large  number  of  active 
channels  that  are  grouped  together.  From  the  time  domain  (Figure  4-5(b)),  we 
see  that  undesired  channels  interfere  with  desired  active  channels.  In  this  case, 
simple  amplitude  or  energy  based  readout  will  fail  the  test.  However,  correlation 
between  those  colored  channels  and  desired  channels  will  still  divide  them  into  two 
different  groups.  Pattern  “1”  can  still  be  perfectly  recovered  by  choosing  the  group 
of  channels  that  have  maximum  readout  in  4-5  (b). 

In  most  cases,  output  level  from  the  undesired  channels  can  be  high  enough 
to  interfere  with  the  desired  active  channels.  Energy  and  distance  based  method 
in  this  case  will  have  difficulty  distinguishing  false  activities  (noise)  while  the 
synchrony  readout  still  extracts  the  phase  information  between  channels  efficiently. 
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Figure  4-5:  Readout  from  an  RKII  network,  (a)  Readout  pattern  that  effectively 
improves  the  discrimination  between  desired  channels  and  colored  channels,  (b) 
Network  response  to  colored  input  pattern  “1” 

However,  instantaneous  detection  of  phase  distributions  may  also  fail  because 
the  frequency  of  oscillations  would  also  be  affected  by  external  inputs.  From 
experiments,  oscillation  frequency  of  the  false  channel  is  close  to  but  not  exactly 
the  same  as  that  of  the  true  active  channels  and  the  phase  difference  between  the 
two  is  time  varying.  That  is,  if  timing  is  not  properly  selected,  instantaneous  phase 
information  will  generate  wrong  information.  However,  from  simulations,  long  term 
correlations  between  two  channels  still  work  because  correlation  between  active 
channels  remains  strong.  The  window  length  used  to  track  correlations  depends  on 
specific  problems.  Nevertheless,  extra  memory  that  stores  previous  estimation  of 
average  values  of  output  state  and  state  powers  will  enable  online  update  of  cross 
correlations  and  save  computational  time  needed  for  long  time-series. 

4.3  Discussions 

Synchronized  oscillations  are  believed  to  play  a key  role  in  many  biological 
models,  especially  for  the  purpose  of  information  processing.  Based  on  our  previous 
work  on  analyzing  how  to  control  an  RKII  network  to  achieve  desired  behaviors, 
we  proposed  to  use  synchronization  as  the  computational  primitive  in  an  RKII 
network.  Information  regarding  collective  behaviors  in  such  networks  is  used  to 
design  the  readout.  Network  parameters  (coupling  strengths  among  different  RKII 
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sets)  are  designed  specifically  to  achieve  desired  responses.  Two  applications  are 
presented.  By  utilizing  two  coupled  RKII  sets  as  a logic  computation  unit,  the 
intrinsic  structure  of  the  coupled  sets  provide  us  with  very  flexible  configurations 
to  implement  any  logic  gates  and  the  system  can  realize  universal  sets  of  logic 
operations.  The  potential  capacity  of  this  network  is  very  large,  however,  the 
controllability  of  all  the  attractors  in  the  output  space  is  not  guaranteed  and  needs 
to  be  further  investigated. 

In  the  second  application,  we  demonstrated  how  to  design  RKII  network  as 
an  associative  memory  by  considering  synchronization  of  groups  of  oscillators 
in  the  output  space.  Finding  parameter  values  in  this  case  is  straightforward 
comparing  to  conventional  learning  algorithms.  Synchronization  as  a readout  uses 
correlation  of  states  to  represent  input  patterns  in  the  output  space.  Particularly, 
the  synchrony  readout  performs  well  in  the  presence  of  noisy  patterns  when  energy 
based  methods  may  fail.  There  are  still  many  problems  that  need  to  be  explored. 
Capacity  and  robustness  of  the  network  in  terms  of  stored  patterns  and  signal- 
to-noise  ratio  need  to  be  investigated.  In  the  applications  pursued  so  far,  only 
fixed  point  solutions  and  synchronized  oscillations  have  been  used  for  encoding 
information.  This  strategy  may  severely  limit  the  capacity  of  neural  network  based 
models.  It  would  be  very  interesting  to  develop  schemes  where  desynchronized 
oscillations,  quasiperiodic  oscillations,  and  possibly  chaos,  together  with  fixed 
point  solutions  and  synchronized  oscillations,  can  all  be  actively  used  to  encode 
information. 


CHAPTER  5 

HARDWARE  IMPLEMENTATION  USING  ANALOG  VLSI 
The  olfactory  system  model  is  a continuous  dynamical  system  described 
by  ODEs.  Analog  circuits  are  the  most  natural  way  to  duplicate  continuous 
computations  for  real  time  applications.  In  this  chapter,  the  RKII  set  is  designed 
using  analog  VLSI  that  works  in  the  subthreshold  region.  Subthreshold  region  is 
used  to  take  into  account  size,  power  consumption  and  functionality.  However,  we 
will  also  see  that  this  may  bring  other  problems  such  as  offset  and  interconnection 
problems.  Based  on  the  architecture  of  an  RKII,  a straightforward  design  is 
divided  into  three  major  tasks  including  summing  node,  nonlinear  function  and  the 
linear  dynamics.  Design  issues  are  discussed  individually.  A new  architecture  that 
combines  the  summing  node  and  the  nonlinear  function  is  also  proposed  to  reduce 
the  complexity  and  power  consumption. 

5.1  Summing  Node 

The  summing  node  receives  external  and  feedback  inputs.  The  output  of  this 
block  is  the  weighted  sum  of  all  incoming  signals.  Specifically,  in  an  RKII,  we 
need  to  implement  a unit  gain  for  external  input  and  user  configurable  coupling 
coefficients  Kmg  and  Kgm.  In  [5],  all  these  weights  were  implemented  using  opera- 
tional amplifiers  outside  the  chip.  External  implementation  has  the  advantage  of 
high  precision,  but  it  cannot  be  scaled  well  and  requires  large  power  consumption. 
Here  transconductance  amplifiers  (transamp)  working  in  subthreshold  region  are 
integrated  on  chip  to  realize  current  additions.  The  coupling  coefficients  are  nor- 
mally greater  than  one,  so  follower  aggregation  circuit  cannot  be  used.  Instead,  a 
current  adder  followed  by  a current-to-voltage  converter  is  implemented  (Figure 
5-1).  Each  of  the  transamps  (Gmi,  Gm2, . . . , GmN ) converts  an  input  voltage  to  its 
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Figure  5-1:  Summing  node  with  N inputs. 


corresponding  current  using  different  transconductances.  Summation  of  all  currents 
is  fed  into  another  transamp  GmB  (base  transamp)  with  negative  feedback  and  is 
converted  into  to  a voltage  output.  Therefore,  the  output  voltage  is 

N 


Vout  = 


i=  1 


mB 


(5.1) 


The  transamp  can  be  designed  simply  using  a five  transistor  differential  pair  shown 
in  Figure  5-2.  When  working  in  the  subthreshold  region,  the  circuit  has  a transfer 
function  as 


Io  = h-h 

eKnVin+/VT  _ eKnVvin_/VT 
(>KnVin+/VT  -j-  gK-nVin— /Vt 

= /6tanh  [^-iVin+  - Fm_)^  , (5.2) 

where  Kn  is  the  capacitive  coupling  ratio  from  gate  to  channel  in  NMOS  transistors 
and  VT  = kT/q  « 25.7  mV  is  the  thermal  voltage.  Kn  is  typically  set  to  0.7.  The 
shape  of  the  I-V  curve  is  sigmoidal.  When  the  differential  input  (V„+  — V^n_)  is 
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Figure  5-2:  Differential  transconductance  amplifier 
small,  (5.2)  can  be  approximated  as 


I0  ~ Gm(Vin+  Fin—)) 


where  Gm  = /(,/«„/ (2 Vy).  It  is  clear  that  this  approximated  linear  relationship 
only  works  when  the  differential  input  is  rather  small.  In  our  case,  linearity  in  the 
summation  block  is  very  important  because  weighting  of  every  input  is  expected  to 
be  a constant.  However,  the  input-output  relations  make  it  hard  to  keep  constant 
weight  values  in  a relatively  large  range  of  input  values.  There  are  many  methods 
to  improve  linearity  such  as  source  degeneration,  bump  circuit  and  floating  gate. 
Figure  5-3  depicts  a transamp  with  source  degeneration  using  double  diffusors  [46]. 
The  output  value  of  IQ  is 


KjT 


2^7 ~(Vin+  - Fin_)  - tanh 


(Vin+  ~ Vin_)  > 

l 1 


I0  — lb  tanh 


2 m + 1 
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Vdd 


Figure  5-3:  Wide  linear  range  tranconductance  amplifier  with  source  degeneration 
using  double  diffusors 

where  m is  the  ratio  between  the  value  of  W/L  of  the  diffusors  (M5-M8)  and  that 
of  the  input  transistors  (Ml  and  M2).  Assume  that  VT  = 25.7mV,  Kn  = 0.7, 
if  linear  range  is  defined  by  1%  distortion  from  the  maximum  transconductance 
obtained  at  Vin+  — Vin _ = 0,  then  the  circuit  shown  in  Figure  5-3  has  an  approx- 
imate linear  range  of  30  mV  [46].  This  is  not  a promising  value  that  could  make 
the  circuit  accommodating  the  large  amount  of  inputs  from  a network.  Based  on 
the  functionality  of  the  transamp,  assume  that  IB  is  the  bias  current  of  the  base 
transamp  GmB,  and  VinjMAX  is  the  maximum  value  of  the  input  signal  before  the 
input  transamp  enters  saturation,  we  have 

N „ 

VinMAx  V 77^  oc  IB.  (5.3) 

~7  ^rnB 
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Moreover,  to  obtain  coupling  strengths  greater  than  1,  Gmi  must  be  greater  than 
GmB  (Figure  5-4).  A large  value  of  output  from  the  series  of  current  adder  could 


Input  Voltage  (v) 

Figure  5-4:  Because  the  coupling  strength  is  always  greater  than  1,  input  current 
range  can  easily  exceed  the  limit  of  the  base  transamp 

easily  exceed  the  linear  range  of  the  base  transamp  in  the  current  to  voltage  con- 
verter, which  makes  the  base  transamp  working  out  of  range  and  the  summation 
inaccurate.  Hence,  the  largest  weight  value  and  input  range  in  this  circuit  are  lim- 
ited. One  way  to  solve  this  problem  is  to  utilize  the  theoretical  analysis  obtained  in 
previous  chapters.  Because  it  is  the  value  of  the  product  of  internal  couplings  that 
determines  the  dynamic  behavior  of  an  RKII,  we  can  average  down  the  individual 
values  of  Kmg  and  Kgm  but  still  keep  the  product  satisfying  the  necessary  operating 
conditions.  Later  in  this  chapter,  a new  design  will  be  presented  that  combines 
the  summation  block  and  the  nonlinearity  so  that  a wider  dynamic  range  and  less 
dependance  on  linearity  can  be  achieved. 
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5.2  Linear  Dynamics 

The  2"d-order  linear  time-invariant  system  can  implemented  using  the  Filter 
and  Hold  (FH)  techinque  [5].  Figure  5-5  shows  the  diagram  of  this  implementation. 
The  two  stages  of  Gm  — C filters  realize  the  functionality  of  a 2nd  order  low-pass 


Figure  5-5:  The  2nd  order  linear  filter.  Filter  and  hold  technique  is  used  to  increase 
the  time  constant. 

filter.  Switches  5j  and  S2  may  suffer  from  charge  injection  and  clock  feedthrough. 
Dummy  switches  Su  and  S2d  are  used  to  neutralize  the  charges  that  may  be 
injected  into  the  capacitor  when  clock  switches  to  ‘off’.  W/L' s of  both  S id  and  S2d 
are  one  half  of  those  of  Sj  and  S2  because  every  time  when  the  clock  switches  to 
‘off’,  half  of  the  channel  charges  from  Si  and  S2  would  be  injected  to  C's. 

Analog  output  of  each  stage  in  the  2nd-order  system  is  held  for  a certain 
period  by  using  the  switches.  When  using  proper  timings  of  the  clocks  $1  and 
$2,  the  actual  time  constant  of  the  system  increases.  The  transfer  function  can  be 
shown  as  [5] 

1 _ e-(Gmi/Ci)kT  j _ e-(Gm2/C2)kT 

= 1 _ e-(Gml/Cx)kTz-\  X 1 _ e-(Gm2/C2)kTz- 1 X Z 

where  k is  the  pulse  width  of  the  clock  (or  inverse  of  the  duty  cycle).  Thus, 
by  setting  proper  duty  cycles  of  the  clocks,  capacitance  (time  constant)  can  be 
effectively  increased.  With  a duty  cycle  of  1%,  the  integrated  capacitors  can  be 
scaled  by  100,  which  greatly  reduces  the  chip  area  needed. 
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The  purpose  of  using  a very  large  time  constant  is  to  mimic  the  exact  fre- 
quencies of  the  measured  data  from  animal  experiments  [1].  However,  in  practice, 
if  all  parameters  are  adjusted  accordingly,  a scaled  time  constant  will  only  affect 
the  speed  of  temporal  development  in  the  network.  Recall  that,  boundary  of  the 
interesting  dynamics  in  an  RKII  is  solely  determined  by  a nonlinear  function  that 
involves  a,  b,  Kmg  and  Kgm.  Therefore,  as  long  as  the  bifurcation  condition  is 
satisfied,  the  overall  dynamics  of  the  network  is  expected  to  be  unchanged.  Actu- 
ally, an  equivalent  scaling  for  both  a and  b does  not  change  the  boundary  set  by 
(a  + b)2/(ab). 

5.3  Nonlinear  Function 

Recall  that  the  nonlinear  function  is  given  as 

{_  ex  — 1 

Qm(  1 ~ e Qrn  ) X > x0  (5.4a) 

— 1 else.  (5.4b) 

In  [5],  the  sigmoidal  shape  is  well  approximated  by  a differential  pair  working  in 
subthreshold  region.  The  asymmetric  shape  with  different  positive  and  negative 
saturations  by  setting  different  sizes  for  transistors  of  the  current  mirror  load  in  the 
transamp.  In  this  way,  Qm  is  predefined  and  can  have  a relatively  high  precision 
if  the  differences  in  transistor  sizes  are  properly  designed  (for  example,  by  using 
parallel  transistors).  However,  ability  to  control  the  actual  shape  of  the  nonlinear 
function  is  also  appealing  to  bring  more  freedom  to  generate  the  desired  dynamics. 
Especially,  different  Q functions  result  in  different  values  of  couplings  necessary 
for  oscillations.  Therefore,  instead  of  building  into  the  chip  a fixed  asymmetric 
ratio,  we  use  an  externally  controlled  current  multiplier  to  adjust  the  current  going 
through  the  load  of  the  differential  pair.  The  circuit  is  shown  in  Figure  5-6. 

We  have 


(5.5) 
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Vdd 


Figure  5-6:  The  nonlinear  block. 

where  I\  and  /2  are  controlled  by  and  Vj,2  respectively.  The  relation  between  the 
differential  input  and  output  current  can  be  expressed  by: 

Ileiyin-VreS)IVT  _ l 

««•)  = frV-'ww  + i (56) 

Equation  5.6  is  similar  to  the  implementation  given  in  [5].  They  both  approximate 
the  sigmoidal  shape  required  by  the  nonlinear  block  in  a KO.  However,  here  the 
actual  asymmetric  ratio  and  shape  of  nonlinearity  in  use  can  be  controlled  by  bias 
voltages  set  outside  the  chip  and  are  not  fixed  at  the  time  of  fabrication. 

An  interesting  observation  is  that  both  the  summation  block  and  nonlinear 
block  use  similar  functional  circuit,  i.e.,  transamps.  The  output  from  the  intermedi- 
ate linear  system  is  shaped  and  added.  Intuitively,  it  seems  feasible  to  combine  the 
procedures  of  shaping  and  summation  to  a single  block.  We  will  show  in  the  next 


section  that  this  is  true  and  can  be  utilized  to  reduce  the  complexity  and  power 
consumption  of  the  whole  system. 
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5.4  Efficient  Design  of  an  RKII  System 

Recall  that  the  general  network  architecture  is  described  by 


1 

ab 


d 2xAt)  . , .d  xAt)  ...  . . 

-^  + (a  + b)-±1  + (ab)xi(t) 


= Ef*  [WvQMt),  Qm ) + Qn),  *)]  + rn 

i,j  = 1,...,N, 


(5.7) 


Current  circuit  design  separates  the  shaping  and  weighting  processes  before  adding 
different  input  and  feedback  signals  together.  In  this  way,  weighting  circuits 
must  satisfy  high  linearity  although  the  signals  being  added  are  all  shaped  by 
a nonlinear  function.  Moreover,  the  current  adder  and  the  nonlinear  function 
are  all  implemented  by  the  same  circuitry,  transamps.  Instead  of  separately 
building  coupling  weights  and  Q(x),  a unified  KmgQ(x)  (KgmQ(x))  function  can 
be  designed  using  transamps  that  realize  weighting  and  nonlinear  shaping  at  the 
same  time.  It  has  the  clear  advantage  of  reducing  system  complexity.  In  addition, 
linearity  is  not  required  anymore  for  the  input  transamps.  The  natural  exponential 
functionality  of  the  differential  pairs  can  be  used  without  any  approximation.  Now, 
the  question  is  how  to  design  a controllable  nonlinear  function  that  can  retain  a 
desired  asymmetric  ratio  and  at  the  same  time  generate  different  power  levels. 

At  least  two  free  parameters  are  needed  to  control  individually  both  the  slope 
and  saturation  level  of  the  nonlinear  function.  Bias  current  can  adjust  both  values 
but  not  separately.  Figure  5-7  shows  a simple  way  to  control  the  ratio  of  the 
current  mirror  in  subthreshold  region. 


Io 


p{vb2~Vb l)  T 
rr 


(5.8) 
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Figure  5-7:  Asymmetric  current  mirror. 


If  vb2  is  set  to  the  power  supply,  output  current  only  depends  on  vbi  ■ Consider  the 
circuit  shown  in  Figure  5-8,  Output  current  can  be  expressed  as 

Ne(Vin+-Vin-)/vT  _ i 


Q(Vm) 


e(Vin+-Vin-)/VT  + 1 ' 


(5.9) 


The  value  of  N is  determined  by  the  bias  Vbr ■ This  nonlinear  function  has  a 


Vdd 


Figure  5-8:  New  implementation  of  the  asymmetric  nonlinear  function 


built-in  offset  at  zero  input.  To  compensate  the  offset,  two  such  transamps  can  be 
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used  [5]  (Figure  5-9).  We  should  note  that  this  scheme  can  only  cancel  systematic 
offset.  Bias  current  determines  the  negative  saturation  where  vbI  can  be  used  to 


control  the  asymmetric  ratio.  Without  the  need  for  shaping  the  output  at  the  end 
of  the  system,  RKII  can  be  redesigned  by  removing  the  nonlinear  function  and 
loosing  the  requirement  of  the  linearity  of  the  adder  (Figure  5-10). 


Figure  5-10:  New  implementation  of  the  RKII  system  that  only  includes  an  adder 
and  a linear  system 

5.5  Digital  Simulation  and  Chip  Measurement 

One  rather  important  intermediate  medium  between  theory  and  analog 
implementation  is  simulation.  In  particular  when  one  is  dealing  with  a distributed 
nonlinear  dynamical  system  for  which  our  analytical  ability  is  still  rather  limited, 
simulations  play  a center  role  because  they  allow  us  to  verify  the  theoretical 
analysis  and  check  the  behavior  of  the  analog  VLSI  chips  we  build. 
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An  equivalent  digital  model  of  the  olfactory  model  is  developed  [5, 47]  using  a 
digital-signal-processing  (DSP)  approach.  The  linear  dynamics  of  the  basic  system 
equation  of  a KO  set  is  described  as 

^x(t)  + ^^x(t)  + x(t)  = P(t).  (5.10) 

The  olfactory  system  is  built  from  this  basic  linear  system  that  is  followed  by 
a static  nonlinearity.  DSP  technique  can  be  used  to  transfer  this  system  to  a 
difference  equation.  The  transfer  function  of  a KO  set  is 


H(s)  = 


ab 


, ,s  + a s + b. 
= ab(- + r). 


(s  + a) (s  + b)  ~Kb  — a ' a — b 
The  time-domain  solution  can  be  expressed  as 

0—bt 


(5.11) 


-,—at 


h(t)  = ab( 


+ 


b — a a — b 


)■ 


(5.12) 


Using  the  impulse  invariant  transformation  technique  with  sampling  period  Ts , the 
discrete  approximation  filter  has  the  impulse  response 


p—anT  p—bnT 

h[n\  = Tsh(nTs)  = Tsab(- 1 -)  (5.13) 

b — a a — b 

Denote  aq  = e~aT,  = e~bT,  C\  — Tsab/(b  — a),  and  ci  = Tsab/(a  — b),  the  transfer 
function  of  the  equivalent  digital  system  is  obtained  as 


H(z) 


Cl  C2 

1 — aqz-1  1 — a2z~l 


(5.14) 


Generally,  the  discretization  of  a continuous-time  system  using  the  impulse- 
invariant  transformation  is  subject  to  aliasing.  However,  the  poles  of  a KO  set 
are  all  positive  and  the  system  has  typical  behavior  of  a low-pass  filter.  Hence, 
aliasing  can  be  effectively  reduced  by  choosing  small  sampling  period  Ts.  Practical 
implementation  of  the  discretized  system  can  be  achieved  using  conventional  ar- 
chitectures such  as  parallel  realization.  When  modeling  a KO  set  specifically  based 
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Figure  5-11:  Comparison  between  impulse-invariant  filter  and  circuit  simulation 

on  the  actual  nonlinearity  designed  in  the  circuit,  we  can  use  digital  simulations 
to  verify  the  circuit  design.  In  Figure  5-11,  bifurcation  boundary  simulated  from 
impulse-invariant  system  is  compared  to  that  obtained  from  circuit  simulation.  The 
result  indicates  a good  match  between  digital  simulations  and  circuit  design. 

The  analog  chip  (Figure  5-12)  was  designed  using  AMI  0.5  technology  and 
includes  an  RKII  set.  Each  KO  set  has  an  area  of  400 x 150 /jm  and  a power 
consumption  of  20  fiW  excluding  the  buffers.  Individual  components  were  also 
fabricated  for  testing  purpose. 

Figure  5-13  shows  the  measurement  of  nonlinear  functions  for  excitatory  and 
inhibitory  cells.  Negative  saturation  is  fixed  at  -100  mV.  Three  different  values 
of  Qm  are  determined  by  the  bias  that  controls  the  ratio  between  positive  and 
negative  saturations.  The  measured  nonlinear  function  has  a relatively  small 
offset  around  5 mV.  Figure  5-14  shows  the  response  from  the  excitatory  cell 
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Figure  5-12:  Chip  layout  of  a single  RKII  set 
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Figure  5-13:  Measured  excitatory  and  inhibitory  nonlinear  functions 

and  inhibitory  cell.  Figure  5-15  shows  the  phase  plot  of  the  measurement.  A 
square  wave  is  connected  to  the  input  of  the  excitatory  cell.  As  we  expected,  the 
qualitative  behaviors  of  the  KII  set  follow  the  theoretical  conclusions  and  exhibits 
two  states  of  dynamical  behaviors  controlled  by  an  external  stimulus. 

As  discussed  in  Chapter  2,  effective  input  is  limited  to  bring  oscillatory 
behaviors.  Taking  a triangular  signal  as  the  input,  the  measurement  (Figure  5-16) 
verifies  that  not  all  positive  inputs  can  induce  oscillation  in  an  RKII  set.  Note  that, 
in  Figure  5-16,  amplitude  of  oscillation  at  a specific  time  is  not  the  steady-state 
amplitude  corresponds  to  its  input  levels  due  to  the  speed  of  the  triangular  input. 
Nevertheless,  the  result  still  demonstrates  the  limitation  on  the  external  input. 

5.6  Discussions 

We  have  designed  and  fabricated  all  the  building  blocks  necessary  for  building 
the  KII  set.  The  dynamical  behavior  measured  from  the  chip  is  qualitatively 
the  same  as  that  obtained  from  analysis  and  simulations.  A chip  includes  8 
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Figure  5-14:  Transient  behavior  measured  from  the  chip.  Channel  1 is  excitatory 
output  and  Channel  3 is  inhibitory  output 

interconnected  RKII  sets  was  also  designed  and  fabricated  (Figure  5-17)  to 
experiment  possible  implementations  when  building  an  RKII  network.  Because 
the  number  of  output  pins  is  limited  to  40,  this  circuit  consists  of  an  analog 
part  (8  RKII  sets)  and  a digital  part  (a  multiplexer).  Despite  the  huge  number 
of  coupling  weights  (bias)  that  need  to  be  set,  only  four  different  bias  voltages 
(weights  obtained  from  outer  product  rule  in  associative  memories)  are  provided 
externally.  The  purpose  of  the  digital  circuit  is  to  multiplex  the  reference  voltage 
so  that  multiple  weighting  circuit  can  share  the  same  one.  However,  the  chip  was 
not  tested  successfully.  In  the  designed  chips,  no  dedicated  circuits  are  used  to 
cancel  offset  and  to  provide  on-chip  bias.  Offsets  of  a fabricated  chip  are  normally 
generated  from  sources  such  as  systematic  mismatch,  fabrication  mismatch  and 
chip  package.  In  a 0.6/im  technology,  a lmV  mismatch  of  the  threshold  voltages 
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is  approximately  equivalent  to  a 300fF  channel  capacitance  at  the  input  [48].  In  a 
highly  interconnected  network,  offset  will  affect  the  performance  in  terms  of  speed 
and  power  consumption.  For  our  purpose,  a deviation  from  mathematical  model 
of  a functional  block  may  directly  affect  the  behavior  of  the  dynamical  system. 

In  a network  of  KII  sets,  where  rich  structure  is  presented  and  probably  with 
sophisticated  boundaries  between  different  attractors,  the  offsets  may  change  the 
desired  behaviors  completely.  To  reduce  the  negative  effect,  electronic  cancelation 
of  offset  should  be  considered  in  the  design.  The  same  idea  applies  to  on-chip 
bias  circuit.  Every  part  of  an  oscillator  needs  to  be  biased.  Most  of  them,  such  as 
linear  dynamics  and  nonlinear  functions,  are  supposed  to  have  exactly  the  same 
functionality  even  in  different  channels.  So  a biasing  circuit  may  be  incorporated 
into  the  chip  so  that  accurate  and  stable  bias  voltage  or  current  could  be  provided. 
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After  these  considerations,  we  should  expect  the  circuit  to  perform  better  on 
matching  the  behaviors  as  what  the  mathematical  models  and  digital  simulators 
can  achieve.  Moreover,  for  every  circuit  simulation  or  observed  offsets  in  the  chip, 
we  could  also  provide  a corresponding  mathematical  model  to  be  used  in  the 
digital  simulation.  In  this  way,  we  have  the  advantages  of  obtaining  fast  simulation 
results  and  verifying  any  theoretical  conclusions  more  efficiently.  This  could  greatly 
reduce  the  effort  needed  to  perform  trial-and-error  procedures  in  the  analog  VLSI 
design  and  increase  productivity.  We  should  also  point  out  that  based  on  previous 
analysis,  qualitative  behaviors  of  a RKII  set  are  mostly  determined  by  the  product 
of  the  two  coupling  coefficients.  Therefore,  qualitative  behaviors  in  the  circuit  of  a 
single  RKII  set  has  relatively  high  tolerance  for  variations  because  the  adjustment 
of  the  couplings  is  straightforward. 
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Figure  5-17:  Chip  layout  of  a mixed-signal  design  that  includes  8 interconnected 
RKII  sets 


When  the  ability  of  the  circuit  to  duplicate  the  results  of  theoretical  models 
can  be  guaranteed,  a KII  or  RKII  network  can  be  built  to  perform  real-time 
applications.  While  we  used  large  time  constants  (large  capacitors)  to  mimic 
biological  system,  when  we  actually  utilize  the  hardware,  a higher  speed  in  the 
part  of  linear  dynamics  may  be  more  desirable  in  terms  of  computational  efficiency. 
At  the  same  time,  the  cost  of  area  in  the  chip  for  the  capacitors  can  be  reduced. 
However,  we  should  also  note  that  the  implementation  of  a KII  network  and  the 
full  Kill  model  in  a chip  also  raise  practical  problems.  The  big  difficulty  is  the 
size  of  the  interconnections  and  weighting  circuits  (Figure  5-17)  required  for  full 
connectivity  (interconnection  grows  with  the  number  of  the  processing  elements 
squared).  A multiplexing  scheme  [49]  or  digitized  signal  communications  such  as 
address  event  representation  (AER)  [50-52]  will  be  needed  to  resolve  this  problem. 


CHAPTER  6 
CONCLUSIONS 

This  dissertation  aims  at  exploring  dynamical  behaviors  of  a kind  of  bio- 
logically inspired  oscillatory  neural  network  (Freeman’s  computational  model 
of  olfactory  cortex)  as  well  as  implementing  such  a network  as  an  information 
processer,  based  on  which  we  hope  to  understand  more  of  the  fundamental  concepts 
of  human  learning  mechanism  in  the  view  of  nonconvergent  networks. 

We  have  obtained  analytical  solutions  for  parameter  optimization  in  an 
RKII  set  and  performed  synchronization  analysis  of  coupled  RKII  sets.  Using 
the  analysis  of  basic  building  blocks  in  the  system,  we  constructed  a system  of 
computational  primitives  based  on  synchronization.  By  using  this  property,  the 
network  can  have  the  ability  to  deal  with  certain  temporal-spatial  patterns  that 
includes  both  amplitude  and  phase  information.  Different  approaches  such  as  PCA, 
energy  based  methods  and  feature  extraction  in  the  frequency  domain  have  been 
studied  in  the  literature.  An  important  job  is  to  determine  which  is  the  best  one 
that  fits  into  the  particular  dynamical  properties  of  our  network.  This  is  still  not 
clear  at  the  current  stage.  Synchronization  is  only  one  of  the  many  rich  dynamics 
that  may  exist  in  the  system.  We  should  notice  that  although  it  shows  potentials, 
it  can  greatly  limit  our  view  from  a global  perspective.  At  the  same  time,  however, 
the  complexity  and  heterogeneity  of  the  brain  discourage  a global  approach. 
Therefore,  we  wanted  to  study  this  hierarchical  model  step  by  step  and  hope  to 
proceed  with  these  initial  analysis  so  that  different  levels  of  complex  behaviors  can 
be  well  recognized  and  controlled  one  by  one. 

More  practically,  from  a system  point  of  view,  there  is  still  a large  extent 
of  unexplored  freedom  in  this  network.  Free  parameters  play  the  center  role  to 
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control  system  dynamics.  At  the  same  time,  they  form  the  bridge  between  the 
input  and  output  spaces  for  obtaining  causality.  For  the  applications  pursued 
so  far,  we  have  not  efficiently  utilized  all  their  capabilities.  Inputs  are  not  time 
varying  and  weighting  is  limited  by  input  space  representation  (coupling  strengths 
is  specifically  designed  according  to  given  input-output  pairs).  However,  we  see 
that  all  the  parameters  bring  diversity  in  dynamics  even  based  on  synchronization 
analysis  where  four  different  behaviors  are  observed  (2  bit  information  in  a 2- 
channel  system).  Different  levels  of  inputs  change  immediately  the  synchronization 
boundaries  obtained  previously  so  that  they  can  possibly  bring  up  all  available  rich 
dynamics  in  a single  fixed-weight  network.  Time  varying  inputs  that  traverse  the 
space  and  break  through  different  attractors  will  certainly  provide  more  freedom  in 
building  the  output  projection  space. 

System  capacity  is  maximized  when  the  richness  of  dynamics  is  obtained 
with  minimal  system  complexity.  Certainly,  the  related  fundamental  question  goes 
back  to  the  problems  of  characterizing  available  attractors  and  representation  of 
output  trajectories  (for  example,  detection  of  trajectories  and  statistical  modeling 
of  transitions  between  attractors).  These  being  said,  a meaningful  future  work 
would  investigate  and  characterize  through  simulations  all  possible  dynamics  in  the 
system. 

from  a general  perspective,  with  the  advances  in  understanding  of  real  biolog- 
ical systems,  mathematical  models  or  engineering  applications  should  also  adjust 
accordingly  and  at  the  same  time  try  to  predict  and  help  discover  new  evidence 
from  biology.  That  is,  the  process  of  building  man-made  biologically  realistic  sys- 
tem should  always  be  an  integrative  process,  where  physiological  studies  provide 
biological  evidence  and  practical  implementations  facilitate  further  understanding 
of  real  biological  systems  as  well  as  help  designing  improved  experiment  to  validate 
engineering  predictions. 
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We  have  designed  and  fabricated  all  the  building  blocks  necessary  for  building 
the  KII  set.  The  dynamical  behavior  measured  from  the  chip  is  qualitatively  the 
same  as  the  simulation.  However,  in  fabricated  chips,  there  are  more  deviations 
and  offsets  in  every  system  block  (due  to  fabrication  process,  design  mismatches 
and  testing  environments,  etc).  Our  analysis  and  chip  measurement  indicate 
that  in  terms  of  qualitative  behaviors,  fault  tolerance  in  the  chip  is  very  positive. 
Particularly,  behaviors  of  an  individual  RKII  is  largely  determined  on  the  product 
of  its  internal  couplings.  This  leaves  much  room  for  choosing  and  adjusting  the 
proper  values  of  control  parameters.  However,  in  a higher  level  where  a network 
of  RKIIs  need  to  be  implemented,  consistency  of  the  performance  is  crucial. 
Advanced  techniques  that  reduce  systematic  offsets  and  noise  should  be  applied 
in  the  future  so  that  variations  can  be  limited.  Moreover,  analog  design  has  the 
natural  difficulty  of  scaling  due  to  the  massive  connections  that  need  to  be  built. 

A complete  DSP  design  or  mixed  signal  design  that  uses  digital  protocols  (such  as 
AER)  for  commutation  among  analog  PEs  may  be  a promising  approach. 


APPENDIX  A 

BASIC  TYPES  OF  LOCAL  BIFURCATIONS 
Bifurcation  is  a change  of  qualitative  behavior  of  a dynamical  system.  Here  we 
introduce  some  basic  types  of  local  bifurcations  [53]. 

A.l  Saddle-Node  Bifurcation 
Consider  a one-dimensional  dynamical  system 


x = f(x,  fi),  x e 5?,  fj,  € 3?",  n > 0,  (A.l) 

where  /i  is  a vector  of  the  bifurcation  parameters,  x0  is  an  equilibrium  at  the 
critical  value  /jl0.  That  is, 

f(x0,Ho)  = 0.  (A. 2) 


If  x0  is  nonhyperbolic  (Equation  A. 3),  / has  a nonzero  quadratic  term  (Equation 
A. 4),  and  / is  nondegenerate  (Equation  A. 5)  with  respect  to  //0,  then  the  system  is 
at  a saddle-node  bifurcation. 


dj_ 

dx 


(z0,//0)  = 0 


av 

dx 2 

dJ_ 

dfi 


(xo,  Ao)  ^ 0 

(xoj  Ah))  ~f~  0 


The  saddle-node  bifurcation  has  a canonical  form  as 


(A.3) 

(A.4) 

(A.5) 


x — a + x2. 


(A.6) 


A. 2 Cusp  Bifurcation 

Consider  a one-dimensional  dynamical  system 

x = f(x,  n),  x e 5R,  n G 9T\  n > 1,  (A. 7) 
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where  p is  a vector  of  the  bifurcation  parameters.  x0  is  an  equilibrium  at  the 
critical  value  hq.  That  is, 


f{x0,Ho)  = 0. 

(A.8) 

If  x0  is  nonhyperbolic  (Equation  A. 9),  / has  a nonzero  cubic  term  (Equation  A.  10) 
but  a zero  quadratic  term  (Equation  A. 11),  / is  nondegenerate  (Equation  A. 12) 
with  respect  to  p0,  and  / satisfies  the  transversality  condition  (Equation  A. 12  and 
Equation  A. 13  are  linearly  independent),  then  the  system  is  at  a cusp  bifurcation. 

^(x0,li0)  = 0 

(A.9) 

93  f 

gx 3 (*0,/h>)  7^  0 

(A. 10) 

d2  f 

dx2  (zo,Ao)  =0 

(A.  11) 

^(z0,/i0)  + o 

(A.12) 

d2f 

The  cusp  bifurcation  has  a canonical  form  as 

(A. 13) 

x = a + bx  + cx 3. 

(A-14) 

A. 3 Pitchfork  Bifurcatation 

Consider  a one-dimensional  dynamical  system 

x = f(x,  fi),  x € K,  fi  € 5ft”, 

(A.15) 

where  /jl  is  a vector  of  the  bifurcation  parameters  and  / is  an  odd  function.  x0  is  an 
equilibrium  at  the  critical  value  /i0.  That  is, 


f(xo,Mo)  = o. 


(A. 16) 
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If  x0  is  nonhyperbolic  (Equation  A. 17),  / has  a nonzero  cubic  term  (Equation 
A.  18)  but  a zero  quadratic  term  (Equation  A.  19),  / is  nondegenerate  (Equation 
A. 20)  with  respect  to  fi0,  and  / satisfies  the  transversality  condition  (Equation 
A. 20  and  Equation  A. 21  are  linearly  independent),  then  the  system  is  at  a pitch- 
fork  bifurcation. 


df 

dx 

d*l 

dx 3 

av 

dx 2 

df 

dfjL 

d2f 

d[idx 

The  pitchfork  bifurcation  is  a special  case  of  the  cusp  bifurcation  when  a 
has  a canonical  form  as 


(x0,Ho)  = 0 

(A.17) 

r(x0,/i0)  7 - 0 

(A. 18) 

|-(afo,  A*o)  = 0 

(A.19) 

(zo,  Ah))  7^  0 

(A. 20) 

(^o,  Ah))  # 0 

(A.21) 

0.  It 


x = bx  + cx 3.  (A. 22) 

A. 4 Hopf  Bifurcatation 

Consider  a two-dimensional  dynamical  system 


x = f(x,y,p)  (A. 23) 

V = f{x,y,(JL)  (AM) 


where  x,  y G 5R,  /jl  G 9ft".  If  the  system  has  a pair  of  purely  imaginary  eigenvalues 
at  x0  and  po>  and  it  satisfies  the  transversality  condition,  then  the  system  is  at  an 
Hopf  bifurcation.  The  most  significant  behavior  of  an  Hopf  bifurcation  is  the  birth 
of  oscillatory  attractor  when  control  parameter  is  changed.  A typical  form  using 
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polar  coordinates  can  be  written  as 

r = ar  + or3  (A. 25) 

0 = w + 7 r2  (A. 26) 

Generally,  if  a stable  limit  cycle  attractor  appears  with  increased  value  of  the 
bifurcation  parameter,  the  bifurcation  is  called  a supercritical  bifurcation.  If 
a the  system  has  an  unstable  limit  cycle  with  decreased  value  of  the  critical 
parameter,  the  bifurcation  is  called  a subcritical  bifurcation.  Detailed  analysis  of 
Hopf  bifurcation  can  be  found  in  [20]. 


APPENDIX  B 

DERIVATION  OF  SUFFICIENT  CONDITIONS  ON  SYNCHRONIZATION  OF 

COUPLED  RKII  SETS 

Eigenvalues  of  matrix  A can  be  computed  analytically  by  solving 


|AI  — A|  = 


A -1  0 0 

ab(l  + Kmm)  A + (a  + b)  — abKgm(Q  (g ^))  0 

0 0 A -1 


—abKmg(Q\m^))  0 ab(l  + Kgg)  A + (a  + b) 


(B.l) 


We  have 


det  ( | AI  — A|)  — [A(A  + a + b)  + ab{  1 + Amm)] 

•[A(A  + a + b)  + ab[  1 + Kgg)] 
-(ab)2KmgKgm(Q\mU))(Q\gM))  (B.2) 

A - ~(a2+  b)  ± [(a  - b)2  - 2 ab(Kmm  + Kgg)  (B.3) 

±2 ab{{Kmm  - Kgg)2  + 4 KmgKgm  ■ (Q' (m^))(Q'(g^))) 1/2]1/2/2  (B.4) 


There  are  two  cases  regarding  the  sign  of  5R[Ai]: 

1.  (Kmm-  I<gg)2  + 4:KmgKgrn(Q'  (g^))  > 0 In  this  case,  the  term  under 

the  outer  square  root  is  guaranteed  to  be  either  real  or  pure  imaginary.  Let 

(a  - 6)2 


Kai  = 


Kb  = 


ab 


KmgKgmiQ'  (m(1)))  (Q'  (tf(1))) 


(B.5) 
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Equations  B.6  and  B.7  are  the  possible  cases  in  which  A has  all  negative 
eigenvalues. 
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(a  - b)2  - 2 ab(Kmm  + Kgg ) + 2 ab^J (Kmm  - Kggf  - AKB  < 0 (B.6) 

(a  - b)2  - 2 ab(Kmm  + Kgg)  + 2 abyj ( Kmm  - Kgg )2  - AKB  > 0 

-(a  + b)  4-  y^(a  - b )2  - 2 ab(Kmm  + A99)  + 2a&^/ (A'mm  - -^ss)2  ~ 4A# 

< 0 (B.7) 


By  solving  Equations  B.6  and  B.7,  we  could  obtain  first  part  of  the  results 
given  in  Section  3.2. 

2.  ( Kmm  - Kgg)2  + 4 Km9Kgrn(Q'(m^))(Q'(g^))  < 0 Let 


R = (a  - b)2  - 2 ab(Kmm  + Kgg) 

I = 2abyJ\ {Kmm  - Kggy  + 4KmgKgm(Q'(mW))(Q'(gW))\ 

a — arctan(— ),  (B.8) 

R 

we  need 

-{a  + b)±Re{y/R  + jI}  < 0.  (B.9) 

Thus, 


Re2{y/R  + jI} 


V R2  + 72cos2  (a/2) 

Lj 

VR2  + 12  + R 

2 

(a  + b)2. 


< 


(B.10) 
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When  Krnm 


I K 99 1 > 2 


Ka2 
2 ’ 


solving  Equation  B.10,  we  obtain 


| (Kmm  - Kgg)2  + 4KmgKgm(Q'(mW))(Q'(gW))\ 
< ~a"ab  + Kmm  + Kgg). 


(B.  11) 


Let 


Ka2 

Kb 


(a  + b )2 
ab 

KmgKm{Q'(mW)){Q’(g^))  , 


(B.12) 


we  have 

(Kmm  — Kgg)2  < AKb  + 2A'j42(2  + Kmm  + Kgg) 

{Kmm  — Kgg)2  > AKB  — 2/Oi2(2  + Kmm  + Kgg).  (B.13) 


From  Equation  B.ll,  we  have 


Kmm  - \Kgg\  > -2,  (B.14) 

so  the  upper  boundary  of  ( Kmm  + |A^99|)2  in  Equation  (B.13)  is  automatically 
satisfied. 

Thus,  we  obtain  the  second  part  of  the  results  in  Section  3.2  as 

0 < Kmm  + |Ap<?|  <--  \J AKb 
-2  < Kmm  ~ \Kgg\ 

{Kmm  — Kgg)2  > AKb  ~ 4LOi2  — 2^42(2  + Kmm  + Kgg).  (B.15) 


APPENDIX  C 

LINEARIZATION  OF  NONLINEAR  FUNCTION 


Simple  nonlinear  functions  can  be  used  to  achieve  similar  input  controlled 


oscillations  in  an  RKII  set.  With  these  modified  nonlinear  functions,  it  may  help  to 
build  digital  circuits  where  only  finite  sections  of  operations  need  to  be  constructed. 
Here  are  two  examples  [54]  of  nonlinear  functions  that  only  have  a few  sections  of 
linear  functionalities. 


One  of  the  key  elements  to  enable  input  controlled  oscillations  in  an  RKII  set 
is  the  maximum  slope  (Q'(x))  of  the  nonlinear  function.  To  achieve  the  desired 
dynamics,  the  nonlinear  function  should  have  at  least  two  different  slopes.  The 
simplest  linearization  is  a three-section  nonlinear  function  as  shown  in  Figure 
C-l.  The  two  slopes  are  Ks  in  the  linear  section  and  0 in  the  saturation  regions. 
Negative  and  large  positive  inputs  result  in  an  infinite  bifurcation  boundary.  Pth 
sets  the  lower  threshold  of  the  input  to  induce  oscillation.  The  upper  bound  of  P is 
determined  by 


where  Qm  is  the  ratio  between  positive  and  negative  saturations.  The  equilibrium 
of  this  system  is  at 


C.l  Example  I 


(C.l) 


m 


9 


* 


P+\KmgKgm\K2Pth 
1 + \KmgKgm\K] 
Kmg  • KS{P  - Pth) 

1 + \KmgKgm\Kl  ' 


(C.2) 


101 


102 


Figure  C-l:  Nullclines  obtained  with  three-section  nonlinear  function 

C.2  Example  II 

Note  that  one  of  the  nullclines  in  the  previous  example  needs  to  have  a lower 
saturation  at  zero  so  that  stable  fixed-point  can  be  obtained  with  zero  input.  The 
approximate  shape  of  the  original  nonlinear  function  can  be  preserved  by  adding 
another  section  in  the  linearized  nonlinear  function  (Figure  C-2).  The  modified 
nonlinear  function  Q(x)  is  defined  by 


Q(x) 


-1 

X < Xi 

(C.3a) 

Ksl  ■ x 

X\  < X < X2 

(C.3b) 

Ks  2 • X 

X2  < X < X$ 

(C.3c) 

Qm 

X < Xs 

(C.3d) 

To  achieve  input  controlled  oscillations,  k 


| Kmg  Kgm  | < 


sl  must  be  less  than  Ks2 ■ When 
1 (a  + b )2 

W,  ’ 


(C.4) 
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Figure  C-2:  Four-section  nonlinear  function 
any  equilibrium  (intersection  between  sections  of  the  tow  nullclines  with  smaller 
slope)  around  zero  is  stable.  When 

i js  „ I . 1 (a  + b )2  ^ ^ 

\KmgKgm\  > ab  , (C.5) 


an  RKII  set  enters  the  oscillatory  state. 
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